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meters.  Only  those  atmospheric  motions  which  are  of  horizontal  scale 
lengths  greater  than  1000  km  and  of  vertical  scale  lengths  greater 
than  10  km  are  considered . 

Annexes  to  this  document  detail  various  factors  to  be  considered 
in  specifying  a  model  of  the  upper  atmosphere. 


Preceding  page  blank 


v 


Abstract 


This  Appendix  to  IDA  Study  S-382  describes  a  numerical  model  for 
the  global  stratosphere,  mesosphere,  and  thermosphere.  The  model  is 
intended  to  be  used  solely  as  a  basis  for  estimating  the  size  of  the 
computational  problem.  The  model  attempts  to  define  the  motions  of  the 
atmosphere,  including  the  troposphere,  up  to  altitudes  of  300  kilo¬ 
meters.  Only  those  atmospheric  motions  which  are  of  horizontal  scale 
lengths  greater  than  1000  km  and  of  vertical  scale  lengths  greater 
than  10  km  are  considered. 

Annexes  to  this  document  detail  various  factors  to  be  considered 
in  specifying  a  model  of  the  upper  atmosphere. 


Preceiling  page  blank 


1  1 

I 

J 


CONTENTS 


Introduction  1 

I.  Coordinates  5 

II.  Kinematic  Relations  7 

III.  Finite  Difference  Mesh  13 

IV.  Variables  17 

V.  Advection  19 

VI.  Horizontal  Divergence  21 

VII.  Eddy  Transport  25 

VIII.  Radiation  Energy  Balance  29 

IX.  Continuity  Equation  35 

X.  Precipitation  39 

XI.  Hydrostatic  Pressure  41 

XII.  Calculation  of  Acceleration  43 

XIII.  General  Plan  of  Calc jiation  45 

Annex  A.  Basic  Conservation  Equations  49 

Annex  B.  Domains  and  Applicability  of  Terms  of  51 

Conservation  Equations 

Annex  C.  Conservation  Equations  in  Spherical  Coordinates  53 

Annex  D.  Conservation  Equations  in  Perturbation  Form  63 

Annex  E.  Difference  Equations  71 

Annex  F.  Considerations  for  Choice  of  Space  and  Time  75 

Steps 

Annex  G.  Heat  Sources  and  Sinks  81 

Annex  H.  Chemical  Reactions  for  Production  and  Loss  83 

Annex  I.  Radiation  Energy  and  Heat  Equations  119 

Bibliography  125 

Preceding  page  blank 

vii 


TABLES 

1.  Even-Odd  Grid  Assignments  18 

2.  Planned  Solar  Energy  Measurements  by  SOLRAD  HI  30 

Satellite  (1973) 

3.  Terrestrial  Radiation  Cooling  Rate  33 

FIGURES 


1.  Horizontal  Scale  Lengths  3 

2.  Spherical  Coordinate  System  6 

3.  Local  Indexing  System  15 

4.  Divergence  Calculation  22 

5.  Atmospheric  Penetration  Depth  of  Energetic  Particles  32 

6.  Calculation  of  Acceleration  43 

C-l.  The  Position  of  the  Rotation  Vector  p  in  the  System  53 

Spherical  Coordinates 

E-l.  Finite  Difference  Approximations  73 

F-l.  Temperature-Altitude  Profiles  of  the  30°,  45°,  60°,  76 

and  75°N.  January  and  Mid-Latitude  Spring/Fall 
Supplement 

F-2.  Components  of  Wind  Velocity  77 

F-3 .  Wind  Near  38°  N  78 

F-4.  Horizontal  Wind  Speed  Spectra  79 

F-5.  Temperature,  Pressure,  Density,  Molecular  Weight  80 


1 


H-l 

Chemical 

Kinetics 

of 

Oxygen  Molecules 

87 

H-2 

Chemical 

Kinetics 

of 

Water  Vapor 

88 

H-3 

Chemical 

Kinetics 

of 

Carbon  Dioxide 

89 

H-4 

Chemical 

Kinetics 

of 

Oxygen  Atoms 

90 

H-5 

Chemical 

Kinetics 

of 

Excited  Oxygen  Atoms 

91 

H-6 

Chemical 

Kinetics 

of 

02  (aV) 

92 

H-  7 

Chemical 

Kinetics 

of 

02  (bXSg) 

33 

H-8 

Chemical 

Kinetics 

of 

Ozone 

94 

H-9 

Chemical 

Kinetics 

of 

Nitrogen  Atoms 

95 

H-10 

Chemical 

Kinetics 

of 

Nitric  Oxide 

96 

H-ll 

Chemical 

Kinetics 

of 

Nitrogen  Dioxide 

97 

H-12 

Chemical 

Kinetics 

of 

Nitrogen  Trioxide 

98 

H-13 

Chemical 

Kinetics 

of 

Nitrous  Oxide 

99 

H-14 

Chemical 

Kinetics 

of 

Hydrogen  Atoms 

100 

H-15 

Chemical 

Kinetics 

of 

Hydrogen  Molecules 

101 

H-16 

Chemical 

Kinetics 

of 

Hydroxl  Radicals 

102 

H-17 

.  Chemical 

Kinetics 

of 

H02  Radicals 

103 

H-l  8 

Chemical 

Kinetics 

of 

Molecular  Nitrogen  Ions 

104 

H-19 

Chemical 

Kinetics 

of 

Atomic  Nitrogen  Ions 

105 

H-20 

Chemical 

Kinetics 

of 

Atomic  Nitrogen  Ions 

106 

H-21 

Chemical 

Kinetics 

of 

Atomic  Oxygen  Ions 

107 

H-22 

Chemical 

Kinetics 

of 

Molecular  Oxygen  Ions 

108 

H-  23 

Chemical 

Kinetics 

of 

Nitric  Oxide  Ions 

109 

H-24 

Chemical 

Kinetics 

of 

Free  Electrons 

110 

H-25 

Chemical 

Kinetics 

of 

Negative  Molecular  Oxygen  Ions 

111 

H-26 

Chemical 

Kinetics 

of 

Negative  Atomic  Oxygen  Ions 

112 

H-27 

Chemical 

Kinet: cs 

of 

Negative  Ozone  Ions 

113 

H-  28 

Chemical 

Kinetics 

of 

Negative  Carbon  Trioxide  Ions 

114 

H-29 

Chemical 

Kinetics 

of 

Negative  Nitrogen  Dioxide  Ions 

115 

H-30 

Chemical 

Kinetics 

of 

Negative  Nitrogen  Trioxide  Ions 

116 

H-31 

Chemical 

Kinetics 

of 

Metal  Atoms,  Metal  Oxide  Molecules  117 

and  Metal  Tons  of 

Na 

,  Mq  and  Cu 

1-1 

1-2 


Equation  Indices 
Indices  by  Domain 


120 

121 


I.  INTRODUCTION 


OBJECTIVE 

The  objective  of  IDA  Study  S-382  is  to  examine  the  feasibility  of 
a  collocation  and/or  consolidation  of  the  Navy  and  Air  Force  current 
and  programmed  (five  years)  operational  environmental  services  com¬ 
puter  centers. 

As  part  of  this  task,  we  present  here  a  numerical  model  for  the 
global  circulation  in  the  strato-,  meso-,  and  thermosphere  (S-M-T), 
showing  features  of  horizontal  scale  larger  than  1000  km.  Since  the 
troposphere  represents  a  .lower  boundary  to  the  upper  atmosphere  des¬ 
cribed,  the  model  describes  the  troposphere  also  in  mesh  coarser  than 
that  commonly  used  for  weather  prediction. 

The  present  model  is  set  up  solely  to  provide  a  basis  for  es¬ 
timating  the  scale  of  the  computational  problem  and  should  not  be  used 
in  any  other  way.  Since  the  described  model,  or  any  other  model  of 
comparable  complexity,  has  never  been  proven  in  practice,  the  details 
may  certainly  be  considered  controversial,  but  such  controversy  need 
not  affect  the  utility  of  the  model  for  the  intended  purpose  of  com¬ 
puter  estimation. 

ORGANIZATION  OF  THIS  DOCUMENT 

The  S-M-T  Model  described  in  the  main  body  of  this  text  includes 
some  but  not  all  of  the  various  considerations  leading  to  a  specifi¬ 
cation  of  a  model  for  purposes  of  its  numerical  computation.  In  Annex 
A,  the  basic  conservation  equations  are  given  in  as  complex'e  a  form  as 
may  be  necessary  to  describe  all  effects,  large  and  small,  that  may 
obtain  in  the  altitude  range  from  sea  level  to  300  km.  In  Annex  B  the 
region  up  to  the  300-km  altitude  is  subdivided  into  several  domains. 
The  terms  described  in  Annex  A  that  are  applicable  to  each  of  these 
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domains,  and  the  boundaries  that  limit  each  domain,  are  also  given. 

In  Annex  C  the  conservation  equations,  described  in  vector  notation  in 
Annex  A,  are  described  in  spherical  coordinates  (r,  0,  \,  t).  In  Annex 
D  the  parameters  described  in  Annex  C  are  given  in  terms  of  the  mean 
quantity  and  a  small  perturbation.  These  relations  may  be  used  to 
treat  the  nonlinear  terms  descriptive  of  the  advection.  In  Annex  C 
schemes  of  difference  equations  are  described.  Leith’s  (1965)  method 
of  treating  advection  is  describeo,  as  is  the  4-point  explicit  scheme, 
with  forward  time  difference,  reviewed  by  Richtmyer  (1367),  which  is 
computationally  stable  if  grid  spacing  and  time  steps  are  suitably 
chosen.  With  computational  stability  and  Nyquist  sampling  consid¬ 
erations  as  given  in  Annex  F,  the  time  and  space  steps  are  defined  to  be 
500  km  horizontally,  5  to  20  km  vertically  and  10  min  in  time.  A 
summary  listing  by  name  of  heat  sources  that  ar-  to  be  considered  in 
each  domain  is  given  in  Annex  G.  In  Annex  H  are  listed  photo  processes 
and  chemical  reactions  described  by  Bortner  and  Kummler  (1968),  as 
significant  in  the  tropo-,  strato-,  meso-,  and  thermosphere.  These 
reactions  lead  to  definition  of  the  production  and  loss  terms  of  the 
number  continuity  equations  and  of  important  heat  contributions  in 
the  atmosphere.  Annex  I  presents  a  timing  analysis  of  radiation 
energy  and  heat  equations. 

OVERVIEW 

Beginning  with  the  effort  of  L.F.  Richardson  (1922),  a  number  of 
attempts  to  compute  the  global  circulation  of  the  atmosphere  have  been 
undertaken.  The  recent  efforts  of  Leith  (1965),  Mintz  and  Arakawa 
(1965),  Smagorinsky  and  co-workers  (1963),  and  Kasahara  and  Washington 
(1967)  have  been  confined  largely  to  simulations  of  the  troposphere, 
problems  of  such  magnitude  as  to  exceed  the  capacity  for  the  fastest 
computers  of  the  day.  The  1971  advent  of  "fourth-”  and  "fifth- 
generation"  computers,  with  capacity  approaching  a  thousand  times 
that  of  computers  used  for  tropospheric  simulation,  promises  that  we 
will  have  the  capacity  to  simulate  the  atmosphere  up  to  and  including 
the  thermosphere,  to  meet  needs  for  data  suggested  by  application  of 
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radio  propagation  for  purposes  of  communication  and  radar,  and  of  orbit 
prediction  of  low-altitude  earth  satellites. 


The  S-M-T  Model  is  one  intended  to  define  the  motions  of  the 
atmosphere,  including  also  the  troposphere  up  to  altitudes  of  about 
300  km,  which  are  of  horizontal  scale  lengths  greater  than  1000  km 
and  of  vertical  scale  lengths  greater  than  10  km. 

In  Fig.  1,  horizontal  scale  lengths  characteristic  of  various 
modes  of  atmospheric  motion  are  described.  The  S-M-T  Model  des¬ 
cribed  herein  will  not  resolve  most  motions  of  turbulence  the  size 
of  hurricanes  or  smaller. 


SECULAR  AND  SYNOPTIC 
PERIODIC  MOTIONS 

GRAVITY  WAVES 

TURBULENCE 

INERTIAL  WAVES 

TIDAL 

FREE  OR 
IMPULSIVE 

GRAV.  TURB. 

INERTIAL  TURB. 

MlCRO-TURB. 

K)\m  lO^km  10km  1km 

3  2 

10  m  10  m  10m  lm 

S5-U-71-*  FIGURE  1 .  Horizontal  Scale  Lengths 

The  method  of  computation  suggested  herein  for  the  S-M-T  Model 
follows  the  pattern  and  has  many  of  the  conventions  described  by 
Leith  (1965).  The  Annexes  detail  various  considerations  leading  to 
a  specification  of  a  model  for  purposes  of  its  numerical  computation: 

The  S-M-T  Model  described  in  the  main  body  of  this  text  includes  some 
but  not  all  of  these  considerations,  neglecting  others  as  elements  of 
negligible  importance. 

In  Annex  A,  the  basic  conservation  equations  are  given  in  as 
complete  a  form  as  may  be  necessary  to  describe  all  effects,  large 
and  small,  that  may  obtain  in  the  altitude  range  from  sea  level  to 
300  km.  In  Annex  B  the  region  up  to  the  300-km  altitude  is  subdivided 
into  several  domains.  The  terms  described  in  Annex  A  that  are  applicable 
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to  each  of  these  domains,  and  the  boundaries  that  limit  each  domain, 
are  also  given.  In  Annex  C  the  bonservation  equations,  described  in 
vector  notation  in  Annex  A,  are  described  in  spherical  coordinates 
(r,  0^  t) .  In  Annex  D  the  parameters  described  in  Annex  C  are 
given  in  terms  of  the  mean  quantity  and  a  small  perturbation.  These 
relations  may  be  used  to  treat  the  nonlinear  terms  descriptive  of 
the  advection.  In  Annex  E  schemes  of  difference  equations  are  des¬ 
cribed.  Leith's  (1965)  method  of  treating  advection  is  described, 
as  is  the  4-point  explicit  scheme,  with  forward  time  difference, 
reviewed  by  Richtmyer  (1967),  which  is  computationally  stable  if  grid 
spacing  and  time  steps  are  suitably  chosen.  With  computational 
stability  and  Nyquist  sampling  considerations  as  given  in  Annex  F, 
the  time  and  space  steps  are  defined  to  be  500  km  horizontally,  5  to 
20  km  vertically  and  10  min  in  time.  A  summary  listing  by  name  of 
heat  sources  that  are  to  be  considered  in  each  domain  is  given  in 
Annex  G.  In  Annex  H  are  listed  photo  processes  and  chemical  re¬ 
actions  described  by  Bortner  and  Kummler  (1968),  as  significant  in 
the  tropo-,  strato-,  meso-,  and  thermosphere.  These  reactions  lead 
to  definition  of  the  production  and  loss  terms  of  the  number  con¬ 
tinuity  equations  and  of  important  heat  contributions  in  the  atmos¬ 
phere.  Annex  I  presents  a  timing  analysis  of  radiation  energy  and 
heat  equations. 
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I.  COORDINATES 


It  is  assumed  in  the  model  that  the  atmosphere  is  always  in 
hydrostatic  equilibrium.  This  means  that  the  pressure  at  a  given 
point  is  determined  by  the  weight  per  unit  area  of  the  air  above 
that  point.  A  differential  expression  of  this  assumption  is  the 
hydrostatic  relation. 


dp  =  -gp  dz  =  -g£mana  dz  (1) 

a 

giving  the  increment  in  pressure  dp  in  terms  of  an  increment  in 
height  dz  for  given  density  p.  Here  g  is  the  (assumed  constant)  ac¬ 
celeration  of  gravity  that  transforms  the  mass  element  p  dz  into  a 
weight  element  gp  dz,  where  p  =X!mana>  the  sum  the  1116155  densities 
of  all  atmospheric  constituents.3 

The  hydrostatic  assumption  rules  out  in  particular  the  possibility 
of  dynamic  pressure  differences  leading  to  vertical  accelerations. 

For  horizontal  scales  of  motion  large  compared  to  the  thickness  of 
the  atmosphere,  this  assumption  is  thought  to  be  quite  valid. 

Another  dependent  variable  that  is  used  is  the  geopotential 
0  =  gz,  giving  the  potential  energy  per  unit  mass  and  serving  as  a 
measure  of  the  height  of  a  given  pressure  surface.  In  the  topography 
of  pressure  surfaces  a  region  of  lower  pressure  corresponds  to  a 
higher  altitude  in  the  z  coordinate. 

The  coordinate  system  by  which  the  atmosphere  is  described  in 
this  report  is  given  in  Fig.  2.  The  earth  is  considered  to  be  a 
sphere  of  radius  a  =  6366  km.  Horizontal  position  coordinates  are 
latitude  0  and  east  longitude  X.  The  vertical  position  coordinate 
is  r  =  a  +  z. 
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II.  KINEMATIC  RELATIONS 


It  is  important  to  distinguish  2  kinds  of  time  derivative.  The 
Eulerian  time  derivative  3/dt  is  based  on  the  time  rate  of  change  at 
a  fixed  point  in  the  space  coordinates  ( 0 »  r),  and  thus  corresponds 

to  the  usual  notion  of  partial 'derivative.  The  Lagrangian  or  sub¬ 
stantive  time  derivation  D/Dt  on  the  other  hand  is  based  on  the  time 
rate  of  change  at  a  point  imbedded  in  ana  moving  with  the  fluid. 

It  is  this  Lagrangian  time  derivative  that  enters  quite  naturally 
in  the  expression  of  many  of  the  physical  laws  governing  the  hydro¬ 
dynamics  of  the  atmosphere. 

The  horizontal  velocity  components  are 

u  —  a  cos  6  ~  (2) 


v  = 


(3) 


positive  toward  the  east  and  pole,  respectively. 


In  the  altitude  coordinate  system  the  vertical  ’’velocity"  com¬ 
ponent  is 


w 


Dr 

Dt 


(4) 


positive  upward.  In  terms  of  these  velocity  components  the  relation 
between  the  two  kinds  of  time  derivative  is 


_D_  _  1 

Dt  dt  u  a  cos  e  3\ 


3  Id.  9 
+  v  —  -r—  +  w  — 
a  30  dr 


(5) 


The  terms  involving  space  derivatives  and  velocity  components  are 
referred  to  as  the  advection  terms. 
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The  equation  of  continuity  is,  for  each  of  the  atmospheric 
constituents  (denoted  by  subscript  a),  as  follows:  . 

iT  *  7  h  <nar2”a)  +  r^Te  &  ‘"aV  *  h  (na  c°s  0  ua> 

=  pa  -  La  (6) 

wherein 

Pa  =  production  of  number  density  of 
d  constituent  (a),  i.e.,  by  dissociation. 

La  =  loss  of  number  density  of  constituent 
a  (_),  i.e.,  by  reconstruction. 

a 

The  conservation  of  mass  requires  the  continuity  relation  of  the 
nondivergent  flow  of  mixture  of  constituents: 

S  +  if  =  0  (7) 


where  &  is  the  horizontal  divergence  of  the  two-dimensional  wind 
field  calculated  in  an  altitude  surface  from  horizontal  components 


u,  v,  i.e., 


a  cos 


1 _  [V 

os  0  [ax 


av  cos  el 


Using  the  boundary  conditions  w  =  0  at  z  =  300  km,  Eq.  7  can  be  in¬ 
tegrated  downward  to  give 


w(z)  - J 


300  km 


£(z)  dz  . 


The  acceleration  equations  in  the  altitude  coordinate  system  are 
the  following,  after  Godske  et  al.  (1957): 
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m 


tan  0  +  —  +  - i -  ^ _  *2 

r  r  cos  0  ax  r  cos  p  a\ 


."a  |  ^rp?  -  TT  -  V  +  ff  +  s  |?  2  n[(sin  u  -  (sin  p  cos  !  )v]  -  sr] 

-  *x[Se  -  ¥ 

-  2  Q^(sin  0  cos  X)  w  +  (cos  p  cos  x)uj  -  S^J 
+  "F"  +  T-  tan  ^  +  7  ^  +  2  o[(sin  X)  w  +  (cos  p  cos  x)v]  -  S^j 


+  e 


=  F  +  F  +  F 
a  em  coll 


(10) 


wherein  Q  =  vector  of  rotation 
0  =  geotential  (gz) 
s  =  specific  volume  =  —  =  1 


p  =  pressure 


p  2>ana 

S  =  force  of  friction  =  V  +  y  c(7  .  v)j 

[4 


=  6 


32(rw)  x  1  a2w  1 


r  39 


a2w 

— 


v  cos  p  ae 


,  tan  0  aw 

+  -7T-ae 


2  3v  2 

ae  +  7? 


r  cos  0 


_3u 

3X 


2w 

“2 


TT  7(v-  V,)] 


+  V 


1  +  5P  [ 


2  tan  0 

- T~  v 


(11) 
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(Landau  &  Lefshitz 
(1959)  15.08,  p.  52) 


The  equation  of  state  of  the  atmospheric  constituents  may  be 
represented  by  the  following: 


P  =  n  kT 
a  a  a 


k  =  Boltzmann’s  constant 

The  thermodynamical  equation  of  energy  conservation  is  written 
as  follows: 


•  •  • 

w  +  s$  =  e  +  ps 


wherein  w  =  rate  of  heat  input  from  sources  radiative  +  conductive 


_  ar  +  7<la''Ta> 


s6  =  rate  of  heat  increase  by  Rayleigh  dissipation  (friction) 


=  unit  volume 

h  /aUi  auk  2  au  \2  I 

“  [2  "  +  5TT  -  3  «ik  +  £<7  *  2) 


=  rate  of  heat  increase  (increase  of  internal  energy 

=  V 


ps  =  compression  energy 


and  symbol  (•)  denotes  +  uy 


In  the  altitude  coordinate  system  the  thermodynamical  equation 
of  energy  becomes: 


10 


/*Ta 
n  C  I  t- - 

i  a  vl  at 


r  cos  9 


+  „  fa  +  ^  +  v  ».) 

Sr  r  cos  e  a\  r  30  / 

'  cl 

>n  ,  n  3n  \ 

r  ^  +  - H - £  +  X  _£) 

3r  r  cos  0  b\  r  36  Ja 

^ ?H(—  SI 


+  w  - —  +  - 

3r  r  cos 


Taw 

h  »SS 

-  1 

A  2  f  3U  \ 

2  ,  Uv  1  sw 

2  3v 

2 

1  2 

aw 

rae 

L 

ar 

T  \3 

r/  7\r  cos  0  3  X/ 

[dr  T30 

?  r3fl 

"  7 

ar 

+  SSL.  +  av  2  an  2  sv_ 

Ir30  r  cos  0  3  X  "  7  r  cos  0  3  x  '  7  rSe 

^  [t^  ^  (rW)  +  r  cos  9  SX  +  r  cos  0  30  ^cos  6  V^J 


wherein  n,  £  are  the  coefficients  of  viscosity. 

The  force  on  a  unit  volume  of  atmosphere  due  to  the  moon’s 
gravity  is  given  in  Jeffreys  (1960),  pp.  231-233. 


/K  M 

r.  I  g  n 


1  3  2 


3.2 


/  j  L  j  ?  C°S  A  W  +  7  Sln  A  Xn  cos  ejJe 

+  [cofl  jl  C0S  A  V  sin  A  \n  (1  +  005  2  e)  |  X 

+  j|  r  sin2  i  xm  sin  2  flV 


wherein  =  X  -  Xm;  Xm  =  longitude  of  moon.  The  electromagnetic 
force  on  a  unit  volume  of  atmosphere  is 
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e  + 
r 


F 

~em 


n  |E 

1  a\°i 


+  E 


u  !a  - 

X  c 


(18) 


The  collisional  force  is 

£colX  =  manavan(^wa  '  V  +  <ua  '  V  \  +  <Va  '  V  S8}  <19) 

where  n  denotes  a  neutral  particle,  e.g.  nitrogen, 
e  is  the  electric  charge  of  species  (,) 

cl  ® 

E  is  the  external  electric  field 
— o 

E  is  the  self-consistent  or  induced  electric  field 
— s 

B  is  the  magnetic  field. 


t 
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III.  FINITE  DIFFERENCE  MESH 


The  finite  difference  equations  are  based  on  2  separate  fixed, 
Eulerian  space-time  meshes  that  are  conveniencly  distinguished  as 
even  and  odd.  Some  dependent  variables  such  as  temperature  and  water 
vapor  concentration  are  defined  on  the  even  mesh  ;  others  such  as 
horizontal  velocity  components  are  defined  on  the  odd  mesh. 

Altitude  z  is  used  as  a  vertical  coordinate;  even  mesh  points 
are  at  41  pressure  levels  (0,  5,  10,  15,  ...,  115,  120  km  in  5-km 
steps;  130,  140,  ...,  180,  200  km  in  10 -km  steps;  and  220,  240,  260, 

280,  and  300  km  in  20-km  steps). 

Latitude  0  and  longitude  \  are  used  as  horizontal  coordinates. 
Even  mesh  points  are  at  even  multiples  of  5°  in  latitude  and  in  longi¬ 
tude  between  the  equator  and  60°  latitude.  Poleward  of  60°  a  coarsen¬ 
ing  of  the  mesh  in  the  longitudinal  direction  is  introduced.  For  65°, 
70°,  and  75°  latitude  the  even  mesh  points  are  at  even  multiples  of 
10°  in  longitude;  for  80°  at  even  multiples  of  20°,  and  for  85°  at 
even  multiples  of  40°.  The  pole  finally  is  an  even  mesh  point.  There 
are  2072  grid  points  in  the  horizontal  mesh  of  a  single  altitude  level. 

The  time  variable  t  is  Greenwich  time.  The  even  mesh  times  are 
even  multiples  of  10  min  starting  from  midnight. 

The  odd  mesh  points  are  at  the  centers  of  the  boxes  defined  by 
the  even  mesh  points.  They  are  thus  at  altitude  levels  of  2.5,  7.5, 
12.5,  etc.,  and  at  times  5  min  later  (or  earlier)  than  an  even  mesh 
time. 

Indexing  of  the  even  mesh  is  by  the  letters  k,  1,  m,  n.  The 
index  k  increases  with  E  longitude  being  0  at  Greenwich  meridian. 

For  latitudes  up  to  60°  the  Greenwich  meridian  is  reached  again  for 
k  =  72,  but  for  example  at  latitude  80°,  it  is  reached  again  for  k  =  18. 
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The  index  1  increases  with  latitude  from  0  at  the  equator  to 
18  at  the  pole.  Note  that  no  distinction  has  been  made  between  N 
latitude  and  S  latitude;  the  two  hemispheres  are  treated  in  the  same 
way,  but  this  implies  that  the  coordinate  system  in  the  southern 
hemisphere  is  the  (left-handed)  reflection  of  that  in  the  northern 
hemisphere. 

The  index  m  increases  (upward)  with  altitude  from  1  for  5  km 
to  40  for  300  km. 

The  index  n  increases  with  time. 

The  odd  mesh  is  then  most  naturally  indexed  with  the  appropriate 
half -integer  indices.  For  purposes  of  problem  organization  and  memory 
storage  assignment,  a  link  has  been  established  between  the  odd  in¬ 
dices  k=^,  1-^,  m-^,  n+^  and  the  even  indices  k,  1,  m,  n. 

The  global  indexing  system  described  is  supplemented  by  a  local 
indexing  system,  shown  in  Fig.  3.  In  this  system,  attention  is  fixed 
on  a  particular  even  mesh  point,  which  is  given  the  single  letter 
index  V  corresponding,  say,  to  the  global  index  (k,  1,  m,  n) .  The 
even  meih  point  neighboring  toward  the  pole  (k,  1  +  1,  m,  n)  is  desig¬ 
nated  ny  P.  The  one  toward  the  equator  (k,  1  -  1,  m,  n)  is  designated 
by  Q.  The  lower  neighbor  corresponding  to  higher  pressure  (k,  1, 
m  -  i,  n)  is  designated  by  U.  The  neighbor  above  corresponding  to 
lower  pressure  (k,  1,  m  +  1,  n)  is  designated  by  T.  The  new  value  at 
the  same  mesh  point  (k,  1,  m,  n  +  1)  is  designated  by  W.  This  leaves 
variation  in  the  longitudinal  direction,  which  is  indicated  by  addition 
and  subtraction;  thus,  the  east  neighbor  (k  +  1,  1,  m,  n)  is  desig¬ 
nated  by  V  +  1,  the  west  neighbor  (k  -  1,  1,  m,  n)  by  V  -  1. 

The  odd  mesh  also  uses  this  local  indexing  system  with  the  same 
letter  designations  through  the  link  described  earlier  between  the  2 
meshes.  Thus  a  westerly  wind  component  u  indexed  by  P  and  written 
as  u^  is  defined  at  the  odd  mesh  point  (k  -  1  +  m  -  n  +  %) . 

The  conceptual  complication  here  is  the  price  paid  for  a  shortening 
of  the  notation. 
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FIGURE 3.  Local  Indexing  System 
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IV.  VARIABLES 


The  independent  variables  whose  subdivision  into  a  discrete 
grid  has  been  described  are  latitude,  9,  longitude,  \,  altitude,  z, 
and  time,  t. 

The  dependent  variables  are  of  2  sorts,  prognostic  and  diag¬ 
nostic.  Of  these  the  prognostic  variables  are  the  more  fundamental 
in  that  they 'define  at  any  given  time  the  state  of  the  atmosphere  and 
its  change.  For  each  of  33  constituents,  there  are  4  variables:  tem¬ 
perature  T,  number  density  n  ,  and  horizontal  wind  velocity  components 
a  toward  the  east  and  v  toward  the  pole,  all  4  being  defined  as  func¬ 
tions  of  latitude,  longitude,  altitude,  and  time.  Of  these  4  variables, 
2,  namely  T  and  n,  are  defined  at  even  grid  points  and  at  even  times; 
whereas  the  other  2,  namely,  the  wind  velocity  components  u  and  v,  are 
defined  at  odd  grid  points  and  at  odd  times.  There  is  1  final  prog¬ 
nostic  variable,  the  surface  pressure,  dp  =  Sn  n  g  dz,  which  is  a 
function  of  latitude,  longitude,  time,  and  altitude.  The  surface 
pressure  pg  is  defined  at  even  latitude  and  longitude  grid  points  and 
at  even  times .  The  prognostic  variables  are  characterized  by  being 
advanced  in  time  through  integration  of  a  first-order  time  differential 
equation  approximated  numerically  by  summation  of  a  first-order  time 
difference  equation.  This  numerical  process  is  of  the  marching  type, 
explicit,  and  involving  only  single  time  differences. 

Any  other  dependent  variables  are  of  the  diagnostic  sort.  They 
are  completely  determined  at  any  given  time  by  the  fields  of  prog¬ 
nostic  variables  at  the  same  time.  The  most  important  of  these  are 
the  horizontal  wind  divergence  &  defined  at  odd  pressure  levels,  even 
latitude  and  longitude  grid  points,  and  odd  times;  the  vertical  "ve¬ 
locity”  component  w  defined  at  even  spatial  grid  points  but  odd  times; 
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and  the  geopotential  0  =  gz  defined  at  odd  pressure  levels,  even 
latitude  and  longitude  grid  points,  and  even  times.  These  assign 
ments  of  variables  to  even  or  odd  grid  points  are  summarized  in 
Table  1. 


TABLE  1.  EVEN -ODD  GRID  ASSIGNMENTS 


I 


Independent  Variable  Grid 

Dependent  Variables 

Altitude 

Latitude 

Longitude 

Time 

Prognostic : 

T 

even 

even 

even 

even 

n 

even 

even 

even 

even 

u 

odd 

odd 

odd 

odd 

V 

odd 

odd 

odd 

odd 

Ps 

even 

even 

even 

even 

Diagnostic : 

£ 

odd 

even 

even 

odd 

w 

even 

even 

even 

odd 

0 

odd 

even 

even 

even 

V.  ADVECTION 


One  of  the  crucial  aspects  of  numerical  models  of  Eulerian  flow 
is  the  approximation  of  advection  terms  in  a  stable  and  accurate  way. 
The  approach  used  in  this  model  is  that  of  Leith  (1965),  which  follows 
the  method  of  fractional  time  steps  used  by  Marchuk  (1964)  and  others 
(Bagrinovsky  and  Godunov  (1957),  D’yakonov  (1962)). 

In  applying  this  scheme  the  difference  equation  can  be  written 

as  a  sequence  of  two  equations,  determining  -|J:  (y>t)  and  -jr  (y,t  +  ^): 

ot  otr 


The  computational  process  involved  in  a  single  time  step  is  divided 
into  2  cycles.  In  the  first  of  these,  advection  in  the  y-direction 
only  is  evaluated  whereas  in  the  second  advection  in  the  x-direction 
only  is  evaluated  starting,  of  course,  with  the  results  of  the  first 
cycle.  The  fractional  time  step  notation  is  a  convenience,  but  the 
results  of  the  first  cycle  have  no  particular  physical  significance. 


VI.  HORIZONTAL  DIVERGENCE 


The  differential  expression  equation  (8)  for  the  horizontal 
divergence  is 

„  _  1  T  9u  ,  3v  cos  0  I 

&  ~  a  cos  '0  [  +  &0  *  J  <8> 

where  a  is  the  (assumed  constant)  radius  of  the  earth. 

In  choosing  an  appropriate  finite  difference  estimate  of  the 
divergence,  it  is  convenient  to  return  to  the  fundamental  definition 
of  divergence  as  net  outflow  per  unit  area.  Figure  4  shows  a  piece 
of  the  grid  in  an  odd  pressure  level  showing  variables  and  indices 
appropriate  to  a  single  divergence  calculation.  The  divergence  &v 
is  that  calculated  for  the  element  of  area  enclosed  by  dotted  lines 
whose  vertices  are  odd  horizontal  grid  points.  Assuming  a  linear 
dependence  of  u,  v  on  position,  the  net  outflow  from  this  element  of 
area  is 


and  the  area  of  the  element  is 


a^sin  81+i 


sin 
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(9) 


Between  2  even  altitude  levels,  winds  and  therefore  also  di¬ 
vergence  are  assumed  to  be  independent  of  pressure,  thus 


w  =  w  ,+$  j(z  -  z  , ) 
m  m-1  m-%  m  m-1' 


with  w-jqq  =  0  in  accordance  with  the  kinematic  boundary  condition  at 
z  =  300. 


VII.  EDDY  TRANSPORT 


The  atmosphere  is  in  fact  in  turbulent  motion.  Any  finite  dif¬ 
ference  description  of  its  state  of  motion  cannot  describe  accurately 
those  scales  of  motion  that  are  small  compared  to  the  grid  size. 

These  scales  are  not  capable  of  explicit  description  and  thus  must 
somehow  be  described  implicitly  or  statistically.  Since  they  are 
not,  then  error  can  result  from  the  natural  transfer  of  energy  from 
the  large  scales  of  motion  to  the  small  scales  of  motion,  which  comes 
about  through  nonlinear  interactions.  As  was  pointed  out  by  Phillips 
(1959),  there  earn  be  aliasing  error  that  in  fact  leads  to  assignment 
of  energy  increments  to  waves  that  have  a  large  rather  than  a  small 
scale  and  thus  falsify  the  long  waves  for  which  the  calculation  pre¬ 
tends  to  be  valid.  One  thing  to  be  done  about  these  small-scale 
eddies  is  to  put  in  some  sort  of  eddy  diffusion  terms  that  would  dif¬ 
fuse  away  singularities  by  damping  short -wave-length  components,  and 
serve  then  as  smoothing  terms.  This  method  of  Leith’s  (1965)  has 

been  used  in  design  of  this  model.  Leith  found  that  a  value  for 

10  2 

the  horizontal  eddy  diffusion  term  coefficient  of  10  cm  /sec  leads 
to  a  stable  calculation.  This  value  is  about  1/4  of  the  gross  Austausch 
coefficient  of  Defant;  i.e.,  should  this  value  have  been  in  the  calcu¬ 
lation  4  times  greater  than  it  was,  presumably  no  interesting  explicit 
eddies  would  have  developed.  The  same  horizontal  eddy  diffusion  co¬ 
efficient  has  been  used  to  diffuse  water  vapor,  temperature,  and  u 
and  v  components  of  the  velocity. 

10  2 

Although  the  value  of  10  cm  /sec  for  the  horizontal  eddy  dif¬ 
fusion  coefficient  was  found  by  experimentation  with  the  Leith  (1965) 
model  itself,  there  is  observational  evidence  that  indicates  that  a 
value  of  about  this  magnitude  is  to  be  expected. 
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L.  F.  Richardson  observed  that  the  eddy  diffusion  coefficient 
appropriate  for  estimating  dispersion  of  clusters  in  the  atmosphere 
was  a  function  of  the  size  of  the  cluster.  Thus,  he  formulated  the 
empirical  law 

D  =  0.2L4/3  cm2/sec  (23) 

2  7 

where  L  is  in  cm  that  hold  over  a  range  of  L  from  10  cm  to  10  cm. 

8  10  2 
For  L  =  2  x  10  cm,  we  would  accordingly  have  D  =  10  cm  /sec. 

According  to  this  empirical  law,  the  largest  scale  disturbances 

in  the  atmosphere  with  wavelengths  of  the  order  of  500  km  would  have 
10  2 

D  =  8  x  10  cm  /sec,  not  much  different  than  Def ant’s  gross  Austausch 
coefficient. 

This  4/3  power  law  of  Richardson  (1926)  has  received  some  theo¬ 
retical  justification  (Batchelor,  1950)  from  the  Kolmogoroff  similarity 
hypothesis  for  equilibrium  turbulent  motions. 

The  introduction  of  a  horizontal  eddy  diffusion  term  in  the 
finite  difference  equations  is  particularly  simple.  The  quadratic 
advection  scheme  (Eq.  20)  used  has  already  introduced  a  second  dif¬ 
ference  term  with  a  dimensionless  coefficient  depending  on  the  ad¬ 
vection  velocity.  To  this  coefficient,  one  need  add  only  another 
dimensionless  term  of  the  form 


„  _  DAt 

r  —  2 
(AX)2 


(24) 


to  include  the  eddy  diffusion  with  the  quadratic  advection  calcula¬ 
tion.  Since  D  is  constant,  r  is  a  precalculated  field  of  values  of 
the  order  of  0.002,  depending  only  on  latitude. 

Also  included  are  vertical  eddy  diffusion  terms  for  most  of 

these  quantities.  The  vertical  eddy  diffusion  coefficient  has  been 

5  2 

chosen  as  2.4  x  10  cm  /sec  at  the  surface  of  the  earth,  decreasing 
as  one  goes  to  higher  levels.  This  value,  determined  more  or  less 
empirically  by  Leith,  produced  an  upward  transport  of  water  vapor  that 

S 
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in  turn  produced  proper  rainfall  and  release  of  latent  heat  in  the 
general  circulation  calculation  made  with  his  model.  The  u  and  v 
components  of  the  wind  are  diffused  vertically  using  the  same  eddy 
diffusion  coefficient  as  for  water  vapor.  However,  entropy  is  not 
diffused  according  to  the  eddy  diffusion  coefficient  prescription, 
for,  if  it  were,  the  heat  flow  would  be  downward,  whereas  it  is  ob¬ 
served  generally  to  be  upward.  This  counter  radiant  flow  of  entropy 
is  of  course  due  to  the  correlation  that  exists  between  the  vertical 
motion  and  excess  temperature.  To  take  this  effect  into  account  at 

t 

least  qualitatively,  a  convection  prescription  has  been  introduced, 
according  to  which  the  vertical  flux  of  heat,  always  upward,  is  made 
a  function  of  the  stability  of  the  atmosphere.  This  convection  function 
decreases  with  increased  stability  for  a  stable  layer,  but  is  constant 

r 

and  equal  to  340  cal/cm*"  day  for  an  unstable  layer.  This  prescription 
permits  a  return  to  stability  in  the  model  atmosphere  without  the 
necessity  of  the  overturning  of  masses  of  air  on  a  scale  comparable 
to  the  mesh  size. 

Neither  the  vertical  eddy  diffusion  nor  the  convection  pre¬ 
scription  permits  eddy  transport  across  the  boundary  defined  by  (p  =  0, 

Z  -  500  km),  which  is  effectively  insulating  and  frictionless.  There 
are,  however,  fluxes  of  momentum,  energy,  and  water  vapor  across  the 
p  =  ps  boundary.  This  boundary  is  at  present  taken  as  a  sea  level  in 
the  model,  with  surface  temperature  Tg  and  the  associated  saturation 
water  vapor  mixing  ratio  |ag  specified  as  a  function  of  latitude  and 
longitude.  Certain  arid  land  areas  are  specified  by  setting  ps  =  0, 
but  no  mountains  should  be  included. 

Similar  flux  formulas  are  used  at  the  surface.  If  we  let 

Vs  =  (u  l  +  v*)*  (25) 

be  the  magnitude  of  the  surface  wind,  and  let  u  and  v  be  the  com- 

s  s 

ponents  of  diffusion  velocity  from  the  surface,  we  can  write  the  flux 
as  follows: 


27 


(26) 


F  =  -C  V  u 
xs  v  s  s 


F  =  -C  V  v 
ys  v  s  s 


where  Cv  is  a  constant  drag  coefficient.  On  the  other  hand,  the  ver¬ 
tical  eddy  diffusion  of  velocity  gives  fluxes  between  the  surface  and 
level  1/2  in  the  amount 


F 


x 


F 


y 


(27) 


The  surface  velocity  components  u  ,  v  are  determined  to  maintain  a 

s  s 

oalance  between  the  stresses  in  Eqs.  26  and  27.  Thus,  we  find 


with 


k 

v 

+  c  V  (p. 


v  s 


(28) 


Similarly,  sensible  heat  flux  from  the  surface  is  calculated 
according  to  a  conduction  formula 

Ft  =  CtVs(Ts  -  T)  (29) 

.-..are  C  is  a  conduccion  coefficient.  Water  vapor  flux  from  the  sur- 

T 

ace  is  calculated  according  to  an  evaporation  formula 


(30) 


where  Cu  is  an 
calculation  of 


C  are 
P 


evaporation  coefficient.  These  fluxes  enter  into  the 
the  changes  in  T  and  p,  respectively.  The  coefficients 
determined  somewhat  empirically  and  are  about  10”'i  gm/cm^. 
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VIII.  RADIATION  ENERGY  BALANCE 


The  heating  of  each  layer  by  incoming  solar  radiation  depends 
on  the  zenith  angle  of  the  sun,  the  energy  absorption  along  the  slant 
path  by  atmospheric  constituents  in  the  layer,  and  the  heat  capacity 
of  the  layer. 

The  very  important  absorption,  scattering,  and  reradiations  of 
solar  energy  by  clouds  and  dust  is  herein  taken  into  account  only 
parametrically ,  in  terms  of  an  effective  and  time-variable  opacity, 
rather  than  explicitly.  This  parametrization  is  particularly  sensi¬ 
tive  to  the  effect  of  clouds  in  the  troposphere  and  may  be  estimated 
on  the  basis  of  cloud  pictures  from  geosynchronous  satellites .  Such 
an  estimation  method  is  not  developed  in  the  present  section,  in  which 
the  older  method  of  Leith  (which  does  not  incorporate  explicit  cloud 
data)  is  described. 

The  solar  declination  6  is  the  latitude  of  the  sun,  a  function 
of  the  day  of  the  year.  The  hour  angle  \  is  the  longitude  of  the 
sun,  a  measure  of  the  hour  of  the  day.  At  a  given  latitude  0  and 
longitude  x>  the  zenith  angle  0  is  given  by 

cos  0  =  cos  6  cos  0  cos(\s  -  x)  +  sin  6  sin  0  .  (31) 

The  sun  is  above  the  horizon  if  cos  0  >  0,  on  the  horizon  if 
cos  0=0,  and  below  the  horizon  if  cos  0  <  0.  The  absorption  of 
electromagnetic  radiation  energy  in  the  band  X-^  to  x2  along  a  slant 
path  of  zenith  angle  X  is 

J±(z)  ni(z)  (ergs  cm'3  sec'1)  (32) 

i 
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(33) 


Ji(z)  =  f2  Gi(X)  I0U>  expj^-F(X)  £  0j(X)  NjCzjJ  dX 
*■1  ^ 


(erg  sec-1  (particle))-1 


I  =  Solar  electromagnetic  flux  outside  the  absorbing  atmos¬ 
phere  at  wave  length  X. 

F(X)  =  Chapman  function  denoting  the  effective  thickness  factor  of 
spherical  shell  atmosphere  =-  sec  X. 


Oj(X)  =  Chapman  cross  section,  a  function  of  wavelength,  of  the  j 

atmospheric  constituent. 

.z 


.th 


N  -  (z)  =  /  n  . 

3  J-znn  1,„  3 


300  km 


(z)  dz  =  integral  number  in  column  along  ver¬ 


tical  path  of  atmospheric  constituent  j . 
-X2 


The  flux 


/*2 


(X)  dx  is  determined  by  SOLRAD  HI  Satellite  (to  be 


launched  in  1973)  for  the  bands  listed  in  Table  2. 


The  deposition  of  the  particulate  energy  of  solar  origin,  as 
measured  outside  the  magnetosphere  from  about  20  kev  to  xOO  Mev  is 
more  difficult  to  determine  as  a  function  of  latitude,  longitude, 
altitude,  and  time  than  is  that  of  the  solar  electromagnetic  radia¬ 
tion. 


The  solar  wind  energy  after  the  SOLRAD  HI  determination  passes 
through  the  magnetosheath  where  it  is  thermalized  by  turbulence  and 
is  rendered  isotropic  (Gosling  et  al.,  JGR  72,  101-112,  1967).  With¬ 
in  the  magnetopause,  the  solar  wind  fills  first  the  outer  zone,  with 
a  spectral  distribution  N(E)  =  k  exp  (-E/Eq).  The  outer  zone,  in 
turn,  empties  into  the  inner  zone,  at  about  10  Rg,  presumably  with 
further  thermalizing  and  isotropizing  due  to  collisions,  and  harden¬ 
ing  due  to  magnetic  pinch.  The  latter  increases  the  population  of 
outer  zone  E  >  40  kev  protons  by  a  factor  of  about  10  ^  (Anderson 
et  al.,  JRG  70,  1039,  1965). 
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TABLE  2.  PLANNED  SOLAR  ENERGY  MEASUREMENTS  BY  SOLRAL  LI  SATELLITE,  1973 


ixeccromagnecic 

Radiation 


Particles 


Band  Coverage 


0.5 

-  3A 

> 

0.5  Mev 

1350 

- 

1550A 

1 

-  5A 

> 

1.0  Mev 

1225 

- 

1350 

1 

-  8A 

> 

2.0  Mev 

c 

•r 

>  1216 

8 

-  16A 

> 

5.0  Mev 

o  1080 

p 

- 

1225 

1 

-  20A 

> 

10.0  Mev 

S’  170 

600 

E 

44 

-  60A 

> 

20.0  Mev 

2  44 

- 

60 

4- 

) 

170 

-  600A 

> 

50.0  Mev 

S  16 

20 

1080 

-  1350A 

> 

100.0  Mev 

w  8 

- 

16 

1225 

-  1350A 

0. 

1  -  0.5A  (20- 

-ISOkev) 

5 

- 

8 

1350 

-  1550A 

0. 

1  -  1.5A  (  7- 

-150kev) 

1 

- 

5 

0.5 

- 

1 

0.1 

- 

0.5A 

20 

-  150  Kcv 

500 

-  1  Mev 

1 

-  2  Mev 

2 

-  5  Mev 

5 

-  10  Mev 

10 

-  20  Mev 

20 

-  50  Mev 

50 

-  100  Mev 

The  consequent  spectrum  is  then  assumed  to  propagate  along  the 
lines  of  magnetic  flux  (along  path  s) ,  to  be  absorbed  according  to 
the  law  of  energy  loss  in  an  electronic  plasma  (Jackson,  [13.88], 
p.  451,  1962),  which  may  be  rewritten  as: 

gf  -  (4nNZe2/m)  yn/--12-3-\)  (34) 

v2  \(kT/m)2/ 

where  the  particles  have  a  charge  number  z,  a  mass  m,  and  a  velocity 
v;  while  the  plasma  has  a  density  N,  a  charge  number  Z,  and  a  tem¬ 
perature  T. 

Cooling  by  radiative  loss  of  heat,  by  emission  from  atomic  oxygen, 
molecular  oxygen,  ozone,  carbon  dioxide,  and  water  vapor,  also  af¬ 
fects  the  heat  content  of  the  atmosphere.  A.  C.  Aiken  and  S.  J. 

Bauer  (1968)  have  found  the  auroral  precipitation  in  geomagnetic 
latitudes  60-70°  to  be  shown  in  Fig.  5. 
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Bates  (1956)  has  described  the  radiative  heat  loss  from  atomic 
oxygen,  by  the  reaction 


to  be 


°(3pi)  -  °(3p2) 


+  hv(\  ~  63[_i) 


_  -elAt 

^  e  ^ 

pq0X  “  noA12El  -En ikT  -E  /kT 

W0  +  W.,  e  +  W  e  °' 


2  1 


-1. 68(10) ~18  e_228/^ 


°  1  +  0.6  e"22E/T+  o.2  e' 


6(10)  n  for  120  km  <  a  <  160  km 


n0  80 


ELECTRON  ENERGY,  KeV 

160  240  320  400  480 


PROTON  ENERGY,  MeV 


FIGURE  5.  Atmospheric  Penetration  Depth  of  Energetic  Particles 
(After  A.C.  A»t»n  nnd  S.J.  Rnuer  1968) 
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wherein 


nQ  =  number  density  of  atomic  oxygen 


=  3  |  statistical  weight 


=  Einstein  coefficient  for  0 


8.9(10)~5  . 


Leith  (1965)  describes  cooling  by  outgoing  terrestrial  radia¬ 
tion  by  specifying  a  cooling  rate  as  a  function  of  pressure  alone, 
independent  of  latitude  or  atmospheric  conditions,  as  shown  in 
Table  3 .  He  hopes  eventually  to  replace  'the  table  by  a  calculation 
of  terrestrial  radiation  balance.  Meanwhile,  the  prescribed  cooling 
rate  interferes  with  a  possible  feedback  control  mechanism  in  the 
model. 


TABLE  3. 

TERRESTRIAL  RADIATION 

COOLING  RATE 

Cooling  rate 

m 

£m 

(deg/day) 

1 

100 

0.21 

2 

200 

1.31 

3 

400 

2.49 

4 

600 

1.74 

5 

850 

0.65 

6 

1000 

0.37 

The  treatment  of  the  radiation  and  heat  source  problem,  from  the  stand¬ 
point  of  the  computer  programming,  is  elaborated  in  Annex  I,  which  also 
gives  estimates  of  central  processor  times  required  for  execution  of 
the  computation. 
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IX.  CONTINUITY  EQUATION 


The  consideration  of  the  equation  continuity  for  each  of  the 
important  constituents  of  the  atmosphere  requires  taking  into  ac¬ 
count  the  mechanisms  for  production  and  loss. 


dn 

_ _ a 

dt 


+  v 


(naHa>  =  P 


(37) 


The  principal  production  and  loss  mechanism,  with  greatest  time 
and  space  variability,  is  that  due  to  ionization  and  dissociation 
resulting  from  the  solar  electromagnetic  radiation,  which  depends 
on  length  and  flux. 

The  dissociation  energies  and  ionization  potentials  of  atmos¬ 
pheric  species  of  interest,  as  listed  by  Bortner  and  Kummler  (1968), 
are  given  in  Annex  H.  The  method  for  estimating  the  effect  of  dis¬ 
sociation  is  illustrated  for  the  production  of  atomic  oxygen  by 
radiation  of  the  Schumann -Runge  continuum 


1300  <  X  <  1750  A 


(38) 


for  example,  is  given  as  follows: 


C>2  +  hv  -*  0  +  0  ;  XD  =  1751  A 

(NHsr  =  PC0)  =  2J(°2)  n<°2)  • 


The  i 


th 


constituent  is  dissociated  according  to 


(39) 


Preceding  page  blank 
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Ji(Xi,X2)  =  Jf^  a±i\)  I0(X)  expj^-F(X)  ^  ajdx(X)  N(z)J  (40) 

wherein 

/* 

n^(z)  dz  for  it  constituent 

t-i 

IQ(X)  =  solar  flux  outside  the  atmosphere 
F(X)  =  Chapman  function  of  solar  zenith  angle  X 
cos  X  =  cos  6  cos  9  cos  AX  +  sin  6  sin  q 
6  =  solar 
9  =  latitude  angle 

AX  =  hour  angle  of  radiation  vector  with  respect  to 
longitude  (GHA  +  X). 

Since  the  radiation  which  is  more  energetic  than  the  Schumann - 
Runge  continuum  may  also  cause  ionization,  it  will  produce  ion  pairs 
at  a  rate  proportional  to  the  next-smaller-integer  value  of  (J/IP) . 

Another  important  source  of  production  and  loss  of  chemical  con¬ 
stituents  is  the  chemical  processes,  defined  by  Bortner  and  Kummler 
(1968)  in  their  Figs.  5-1  to  5-31,  inclusive,  which  are  also  described 
in  Annex  H.  Only  the  most  important  reactions  and  their  effects  are 
shown.  The  species  are  arranged  in  the  following  order: 

a.  Major  neutral  species:  H20,  CO^,  0 2 

b.  Minor  neutral  species:  0,  0(1D),  02(1A),  02(1e),  0^, 

N,  NO,  N02,  N20,  H,  H2,  OH,  H02, 

c.  Positive  ions:  N+,  N+,  0+,  0*,  N0+ 

d.  Negatively  charged  species:  e,  0".  0",  0~,  C0~,  NO",  NO" 
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e.  Metallic  species:  of  form  M,  MO,  M02,  M+  of  metals,  Na, 

Mg,  and  Cu. 

For  each  of  the  reactions  shown  in  Appendix  H,  a  rate  coefficient  is 
given  in  centimeter-gram-sec  (cgs)  units,  that  is,  in  cm3  sec  1  per 
particle  for  two  body  and  in  cm  sec"  per  (particle)  for  three  body 
reactions.  Particle  densities  are  in  cm'3  and  temperature  in  °K. 

M  is  used  to  represent  a  collision  partner.  All  values  for  rate  con¬ 
stants  are  given  by  the  constants  a,  b,  and  c  for  a  rate  constant 
of  the  form  R  =  exp(-c/T).  The  rate  coefficients  of  Appendix 

H  may  have  error  factors  of  2  to  10. 

Bortner  and  Xummler  (1968)  list  33  atmospheric  species,  with 
approximately  8  reactions  each,  or  244  reactions.  The  reactions  are 
to  be  accounted  in  each  of  33  equations  of  the  form  of  Eq.  6,  in 
general,  or,  for  example,  following  Fig.  H-4,  of  that  for  n(0): 


^St21  =  2  (mop-  (X<1750A>)  n(°2>  -»  (mo^) (X<3105A))  n(°z 


^£^1)(X<3981A))+  2Ri  n(02+)  n(e)  -  RjN(0)n(02)n(M)  -  RRn(02)n(M) 


(41) 


-R1n(0)n(0H)  -  RMn(0)n(H02)  -  Rnn(0)n(03) 


where  the  J/IP  expressions  take  the  value  of  the  next  smaller  integer. 
The  33  continuity  equations  are  solved  simultaneously  at  each  grid 
location . 
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PRECIPITATION 


Tentative  new  values  of  temperature,  T,  and  water  vapor  mix¬ 
ing  ratio,  p.  =  nw/^ni  are  calculated  by  Leith  (1965)  ignoring  the 

possibility  of  condensation  and  associated  release  or  latent  heat, 
but  taking  into  account  all  other  sources  and  sinks.  This  tentative 
thermodynamic  state  is  then  examined  to  see  if  it  corresponds  to  a 
state  of  supersaturation.  If  not,  these  tentative  values  of  T  and 
U  are  taken  as  the  final  values.  Otherwise,  condensation  occurs 
leading  to  a  decrease  in  p  and  an  increase  in  T  to  establish  satura¬ 
tion. 

Let  S(T)  -  Pus(T)  be  the  tabulated  function  proportional  to 
vapor  pressure  defining  saturation  values  ps  of  p;  then  the  adjust¬ 
ment  calculation  is  based  on  the  following  two  equations 


q'/n  -  SCT11*1)  -  S(P) 
b  u;  "  Tn+1  _  p 

and 

Cp(Tn+1  -  T)  =  L(nn+1  -  p) 

where 

S(T  )  =  py-s(T  )  =  pp. 


(42) 

(43) 


and  L  is  the  latent  heat  of  vaporization. 

We  solve  this  pair  of  equations  in  the  2  unknown  quantities 
T0*1,  qn+1  by  eliminating  Tn+1  -  T  ;  we  have 

-S'(T)-jf-  (Un+1  -  F>  =  P^+1  -  S(F)  (44) 

P 

Preceding  page  blank 
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giving 


Here  L/C 

P 


,  S(T>  +  (L/C  _)S'(W 

n+1  <  _ 

jT  +  (L/Cp)S'(f) 

Tn+1  =?+  L  (~  _  ^n+i)  # 

P 

=  2483  deg. 


(45) 


(46) 


XI.  HYDROSTATIC  PRESSURE 


The  hydrostatic  integration  can  be  carried  out  after  a  new  value 
of  surface  pressure  is  obtained. 

The  downward  integration  of  the  horizontal  divergence  to  de¬ 
termine  uu  terminates  at  the  surface  pressure  pg  with  a  surface  value 

uu_  • 
s 

From  the  kinematic  surface  boundary  condition  as  given  by 
Eq.  9,  we  have  the  surface  pressure  evolution  equation 


The  advection  involves  surface  wind  components.  It  is  evaluated  as 
all  other  terms  by  quadratic  interpolation  in  separate  sweeps  for 
north-south  and  east-west  contribution.  The  pressure  at  altitude 
Z  may  now  be  determined  by  integrating 


P(Z)  =  ps 


(48) 
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XII.  CALCULATION  OF  ACCELERATION 


The  acceleration  term  coming  from  the  gradient  of  pressure  is 
approximated  as  follows: 


(49) 


The  centering  of  these  estimates  can  be  seen  in  Fig.  6. 


<  0.I1.T1.I 

FIGURF  6  Calculation  of  Acceleration 


Preceding  page  blank 
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The  Coriolis  terms  are  evaluated  at  the  new  time;  thus,  if 


F  =  ^2n  sin  e  +  - - At 


(50) 


the  acceleration  equations  become 

un+3/2  .  „n+l/2  .  fi  un+l/2  +  Fo™-3/2  .  M 


a  v 


At 


n+3/2  n+1/2  n+1/2  _  . 

uv  -  uv  =  ^auv  “  Fu 


•  (■»). 


(51) 


At 


where  A  indicates  the  change  due  to  advection  and  diffusion. 

cl 


if 


now 


_  n+3/2  n+1/2 

G  =  uv  +  Aauv 


■¥l 


At 


(52) 


=  un+3/2  +  A  Un+1/2  -  (m 


H  =  %  '  +  vr'  -  (tjv  ^ 


we  have  finally 


n+3/2  _  G  +  FH 

Uv  "  IT? 


(53) 


u 


n+3/2  _  H  -  FG 


1  +  F 
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XIII.  GENERAL  PLAN  OF  CALCULATION 


The  way  in  which  the  various  influences  described  are  assembled 
into  a  single  complex  model .follows  the  method  of  fractional  time 
steps  used  by  Marchuk  (1964)  and  described  in  the  discussion  of  ad- 
vection. 

The  prognostic  variables  are  assumed  to  be  known  at  a  given  time; 
the  calculational  cycle  for  1  time  step  advances  these  variables  to 
new  values  for  a  time  10  min  later.  It  is  assumed  at  the  beginning 
of  the  calculational  time  step  that  the  u,  v,  components  are  known 
at  the  odd  time  1/2  time  interval  or  5  min  later  than  the  time  at 
which  the  variables  n  ,  T,  p  are  known. 

cL  S 

The  calculation  for  a  complete  time  step  is  broken  down  into 
separate  cycles.  The  first  of  these  computes  by  Eq.  20a  new  values 
of  all  dependent  variables  based  solely  on  advection  in  the  north- 
south  direction.  Using  these  results,  the  next  cycle  computes  by 
Eq.  20b  new  values  based  solely  on  advection  in  the  east -west  direc¬ 
tion.  The  coefficients  fcr  the  advection  calculations  for  the  even 
mesh  variables  n  and  T  are  obtained  by  averaging  the  8  surrounding 
values  of  velocities  defined  on  the  odd  mesh.  All  the  rest  of  the 
contributions  to  change  in  the  values  of  the  dependent  variables 
are  computed  in  a  final  cycle. 

In  this  final  cycle,  one  starts  at  the  top  of  the  atmosphere 
and  evaluates  by  Eq.  21  successively  working  downward  the  horizontal 
divergence  &  of  the  wind  at  odd  altitude  levels  but  at  even  hori¬ 
zontal  mesh  points.  The  negative  of  this  divergence,  -£,  is  inte¬ 
grated  downward  from  the  top  to  give  by  Eq.  22  the  vertical  "velocity" 
w  on  the  basis  of  the  continuity  equation  describing  conservation  of 
mass.  In  this  way  the  diagnostic  variable  w  is  defined  at  even  mesh 
points. 
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A  tentative  new  value  of  T  is  obtained  by  satisfying  the  finite 
difference  approximation  to  the  law  of  conservation  of  energy  Eq.  16 
taking  into  account  the  vertical  advection  terms  as  well  as  any  pos¬ 
sible  energy  sources  or  sinks  but  excluding  temporarily  the  source 
coming  from  the  release  of  latent  heat  associated  with  possible  con¬ 
densation.  Similarily  a  tentative  new  value  of  |j  is  obtained  by 
taking  into  account  vertical  advection  of  constituents  in  satisfy¬ 
ing  the  finite  difference  approximation  to  the  law  of  conservation 
of  numbers  (Eqs.  6,  37)  taking  into  account  any  sources  and  sinks 
but  ignoring  for  the  moment  the  possible  loss  of  water  vapor  by  con¬ 
densation  and  precipitation.  Should  the  resulting  tentative  new 

values  of  n  =  n(Ho0)/£n  and  T  represent  an  atmospheric  state  in  a 
^  a  a 

condition  of  supersaturation,  then  using  Eqs.  45  and  46  the  appro¬ 
priate  amount  of  water  vapor  is  removed  and  the  corresponding  amount 
of  energy  is  introduced  to  return  the  atmosphere  to  exact  saturation. 
The  resulting  values  of  p.  and  T  are  taken  as  final  for  that  time 
step.  In  the  other  case,  that  is  of  subsaturation  conditions  after 
the  tentative  u,  T,  calculation,  then  the  tentative  results  are 
taken  as  final.  This  calculaticnal  process,  which  results  in  the  new 
values  of  p  and  T  proceeds  from  layer  to  layer  down  through  the  at¬ 
mosphere  to  the  bottom. 


At  the  bottom  by  carrying  the  integration  of  -8  to  the  limit 
p  the  value  of  w  appropriate  to  the  surface  is  obtained.  This 

o  S 

together  with  an  extrapolated  value  of  the  wind  fields  to  give  an 
estimate  of  the  surface  wind  components  u  ,  v  permits  an  evaluation 

o  o 

of  the  new  surface  pressure  through  the  pressure  tendency  Eq.  47. 

Knowing  the  new  value  of  surface  pressure  and  knowing  from  sur¬ 
face  topography  the  geopotential  0s  to  which  it  relates,  it  is  pos¬ 
sible  to  integrate  the  hydrostatic  Eq.  47  upward  now  making  use  of 
the  newly  acquired  values  of  water  vapor  and  temperature  to  find 
new  values  of  the  pressure  at  the  odd  altitude  levels.  These  new 
values  of  pressure  are  known  1/2-time  step  later  than  the  u,v,  com¬ 
ponents  on  the  same  altitude  levels  but  at  the  surrounding  even  hori¬ 
zontal  mesh  points.  This  permits  their  gradients  to  be  appropriately 
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centered  for  purposes  of  evaluating  the  acceleration  term  in  the 


momentum  Eq.  51.  For  evaluating  the  vertical  advection  terms  in  the 
momentum  equation,  it  is  necessary  to  average  surrounding  values  to 
find  a  w  appropriate  to  the  odd  mesh  point  at  which  u,  v,  components 
are  defined.  In  the  u  momentum  equation,  the  v  entering  into  the 
Coriolis  term  is  taken  at  the  new  time;  and  for  the  v  equation  the 
Coriolis  term  value  of  u  is  taken  at  the  new  time.  This  represents 
an  evaluation  of  the  Coriolis  terms  1/2 -time  step  too  late  and  is 
done  as  a  practical  macter  to  provide  a  small  amount  of  damping  of 
inertial  oscillations  in  the  model.  This  calculational  process  of 
evaluating  new  values  of  p  by  integration  of  the  hydrostatic  Eq.  48 
and  evaluation  of  new  values  u  and  v  by  time  integration  of  the  mo¬ 
mentum  Eq.  51  proceeds  upward  through  the  atmosphere.  Finally,  when 
the  top  of  the  atmosphere  is  reached,  again  1  time  step  has  been  com¬ 
pleted  in  the  calculation;  and  new  u,  v,  Na,  T,  and  p  values  have 

s 

been  determined. 


47 


ANNEX  A 


BASIC  CONSERVATION  EQUATIONS 


mckehu*  (i,  ii,  nij 

Va(s?  +  H,  •  7H.)  =  £ a  -  «V"»  +  W  +  n,«a(5o  +  ^  ' 

(1)  (2)  (3)  (4)  (5) 


“  mana  ^  *  Hb^  +  v^C^2 Qj  +  9  •  u^)  +  raying.*  mjnju^  x  Q 


(6) 

* 
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rate  of  increase  of 
internal  energy 


isr  ♦  s.  ■  J 
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Temporal 

Change 

Advection 

Compression 

Rate  of  heat 
increase  by 
conductivity 

Rate  of  heat  increase 
by  Raleigh  dissipation 

EQUATION  OF  STATE  (VI) 

p  =  nk.T 


rate  of  increase  from 
radiative  and  conductive 
heat  sources 

1$=) 


(6) 

Rate  of  heat  input 
from  heat  sources  and 
radiation 


Preceding  page  blank 


49 


ANNEX  B 

DOMAINS  &  APPLICABILITY  OF  TERMS  OF  CONSERVATION  EQUATIONS 
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ANNEX  C 

CONSERVATION  EQUATIONS  IN  SPHERICAL  COORDINATES 


The  following  description  of  the  conservation  equations,  in  the 
form  of  the  spherical  coordinate  system  employed  for  the  S-M-T  Model, 
is  adapted  from  that  of  Godske  et  al.  (1957),  for  the  mixture  of  all 
constituents  as  one  fluid. 


THE  EULERIAN  EQUATIONS  OF  MOTION  AND  CONTINUITY 

In  spherical  coordinates  (r,  \,  ep)  as  given  in  Fig.  C-l,  sup¬ 
pose  the  vector  of  rotation  Q  to  be  situated  in  the  meridian  plane 
X  =  0,  and  to  form  an  angle  cpQ  with  the  equatorial  plane  cp  =  0. 
Then,  the  components  of  Q  are, 


(C-l) 


S5-1J-71-7 

FIGURE  C-l .  The  Position  of  the  Rotation  Vector  O  in 
the  System  Spherical  Coordinates 

blank 


Preceding  pa?n 
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We  introduce  the  abbreviations 


i 


(®)  ■fe-rfe 


^  =  hi  0-r.  + 


r  cc 

1  S 


5r  ’"r  r  cos  to  SX  X  r  5cp  to 


+  JSI. 
r  Sep 

(C-2) 

l  a 

7  0i 

r  ocp  CD 

(C-3) 

_  1  a 


bv. 


a(v.  cos  cp) 


v.v  =  —x  -2—  ( r  vr)  +  -  — + - -A— 

-  p2  Sr  r  cos  tp  a^  r  cos  tp  Sep 


(C-4) 


vxv  =  — - f 

—  r  cos  co 


5v 

air-  -  ax  (cos  * 


v]  '^p-b  <~v] 


JP 


r  cos 


9 


ai  (r  cos  <p  H 


5v 

"ax 


(C-5) 


0  ~  V- 

r 


cos  tp  [^r(rcos  9  &*)  +  ^  (cos  a\) 

h  (cos  *  If)] 


(C-6) 


and  the  nomenclature 


Q  =  vector  of  earth's  rotation 
v  =  velocity 
0  =  geopotential  (gz) 
s  =  specific  volume 

5  =  force  of  friction 
p  =  pressure 

p  =  density 
9  =  temperature 

6  =  dissipation 

e  =  internal  energy  per  unit  mass 
w  =  heat  per  unit  mass  from  sources 
F  =  lunar  tide  force 


? 
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ml 


3  o (  o  o  2 

=  pv  •  r  -4—  r  3  -  cos  Xm  +  sin  i  cos  cp 


£em  =  electromagnetic  force  ( 

Icon  =  collision  force  =Evab(^  -  vb) 

b 

P-L  =  rates  of  production  minus  loss  of  number  density 


(C-8) 


(C-9) 


we  find  for  the  equation  of  motion 

‘rlar),  -  \  \  -  If  ♦  5 B  - 2n  cos  ■*>  sln  x  v» 


-  2q  (cos  cp  sin  \  -  sin  cp  cos  cpo  cos  X)v^  -  Sr 


p. 


V  V  V  V  „ 

X  cp  ,  r  X  ,  1  .  s  dP 

r  r  r  cos  cp  dX  r  cos  cp  dX 


+  2(1  (sin  cfe  cos  cp  -  cos  cp0  sin  cp  cos  \)rv 

2n  (sin  cpj  sin  cp  +  cos  cft>  cos  cp  cos  \)v^  "  S^J 

2 

f/dv  \  v  v  v  , 

^[(ar2) s  +  -f-1  +  4" tan  *  +  7  U  +  2n  003  sin  x  vr 


C-10) 


+  2Q  (sin  cp-.  sin  cp  +  cos  cpb  cos  cp  cos  \)v  -  S 


=  F  +  F  +  F  ,, 
—  — em  —coll 
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Similarly,  the  equation  of  continuity  can  be  written  in  the  form 

-2  as  13,-12,  1  SS  V> 

at  2  3r  ^  r'  r  cos  cp  a\ 


(cos  mS  v  )  —  P  —  L 
r  cos  cp  3cp  ™  cp 


(C-ll) 


<Lg  +  i_  (pr2  v  )  +  - i -  |_  (pv  )  +  - - - (p  cos  cp  v  )  =  P-L 

at  T  Jr  ^  r  cos  cp  3\  X  r  cos  cpdcp  y  cp 

r  (C-12) 

The  equation  of  motion  is  essentially  simplified  if  n  vanishes,  or 
is  directed  along  the  axis  of  the  system  of  coordinates,  so  that 
cpo  =  tt/2. 


THE  BOUNDARY  CONDITIONS 

The  equations  of  motion  and  continuity  apply  directly  only  to 
the  interior  of  a  fluid.  At  external  boundaries,  or  at  internal 
surfaces  of  discontinuity,  the  solutions  of  those  eau<  tions  have  to 
satisfy  certain  boundary  conditions.  We  shall  here  consider  only 
the  frictionless  case. 


f(x,  y,  z,  t)  =  0 


(C-13) 


denote  the  equation  of  a  surface  of  discontinuity,  generally  mov¬ 
ing  in  space,  but  at  rest  if  t  does  not  appear. 

The  kinematic  boundary  conditions , 


|f  +  v  •  vf  =  0  , 


+  v'  •  vf  =  0  , 


(C-14) 


express  the  fact  that  the  p articles  of  velocities  v  and  v'  on  both 
sides  near  the  surface  of  discontinuity  will  follow  this  surface 
during  its  motion. 

Formula  C-14  can  also  be  written 
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v  •  Vf  =  V7  •  7f  =  -  ■— T  ,  (C-15) 

ot 

an  equation  that  states  that  the  velocity  satisfies  the  solenoidal 
boundary  condition  expressing  the ^continuity  of  the  normal  velocity; 
in  a  sense  it  is  the  form  taken  by  the  equation  of  continuity  near 
a  surface  of  discontinuity. 

The  dynamic  boundary  condition, 

p(x,  y,  z,  t)  =  p'(x,  y,  z,  t)  ,  (C-16) 

expresses  the  continuity  of  the  pressure  at  the  surface  of  discon¬ 
tinuity,  and  may  thus  be  said  to  represent  the  equation  of  motion 
at  such  a  surface. 

In  many  applications  the  mixed  boundary  conditions 

P  )  +  v-v(p  -  p7)  =  0,  — -Pa~  P  ^  +  v'.vCp  -  P7)  =  0  (C-17) 

found  by  combining  the  kinematic  and  dynamic  conditions  are  especially 
useful. 

The  conditions  at  a  free  surface  where  the  pressure  vanishes 
are  obtained  by  simply  annulling  all  primed  letters  in  the  above 
expressions.  The  condition  at  a  rigid  surface  is  represented  by 
the  first  of  equations  C-14  stating  that  the  normal  velocity  of  the 
fluid  is  equal  to  the  known  normal  velocity  of  the  boundary. 

The  above  boundary  conditions  can  easily  be  written  explicitly 
in  a  component  form  in  any  system  of  orthogonal  coordinates. 

THE  EQUATIONS  OF  STATE  AND  ENERGY 

It  can  be  shown  by  physical  experiments  that  the  variables  of 
state — pressure,  temperature,  and  density --are  mutually  related 
by  the  so-called  equation  of  state.  This  equation  can  be  shown  to 
be  independent  of  the  motion  (see,  for  instance,  0.  Bjorgum,  1944, 
pp.  26-27).  The  equation  of  state  may  be  fairly  easily  expressed 
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in  Langrangian  coordinates  that  refer  directly  to  the  moving  particles 
If,  however,  the  Eulerian  coordinates  (x,  y,  z)  are  used,  the  problem 
of  expressing  the  equation  of  state  by  means  of  these  coordinates 
arises.  The  solution  of  this  problem  is  generally  very  difficult, 
but  becomes  simple  for  an  atmosphere  whose  particles  do  not  change 
their  composition. 

Combining  the  equation  of  state  with  the  equations  of  motion 
and  continuity,  we  arrive  at  a  system  of  5  scalar  equations.  In 
general  this  system  is  incomplete,  for  there  are  6  dependent  variables 
It  is  therefore  necessary  in  the  general  case  to  introduce  the  energy 
equation  as  a  sixth  relation. 

EULERIAN  FORM  OF  THE  SIMPLIFIED  EQUATION  OF  STATE 

Ir.  the  Eulerian  system  of  coordinates,  the  values  A,  B,  ...  of 
the  coordinates  of  numeration  are  unknown  until  all  integrations 
have  been  carried  through.  The  Lagrangian  form  %P  -  DR0  =  0  of  the 
equation  of  state  is  therefore  impractical.  An  equation  of  state 
adapted  to  the  Eulerian  method  is  found  through  elimination  of  the 
quantities  A,  B,  ...  This  elimination  is  effected  by  means  of  the 
Lagrangian  equation  and  those  obtained  from  it  by  individual  differ¬ 
entiations,  namely, 

X  =  | £  +  X=(|t  X,...  (C-18) 

In  the  general  case  where  the  parameters  A,  B,  ...  are  variable 

•  •  •  •  •• 

with  time,  new  unknown  variables  A,  B,  ...,  A,  B,  ...  are  thereby 
introduced,  and  it  is  necessary  to  add  more  equations,  namely  those 
which  characterize  the  individual  variation  of  A,  B,  ...,  e.g.,  the 
variation  of  water  vapor.  When,  however.  A,  B,  ...  are  assumed  to 
be  individually  constant,  A,  B,  ...,  A,  B,  ...  all  vanish.  Then 
A,  B,  ...  can  be  directly  eliminated  and  a  Eulerian  equation  of  state 
can  be  obtained  having  the  following  form: 

X(P>  *»,  Q,  p,  s,  9 ,  p,  ...)  —  0  (C-19) 
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The  equation  of  state  is  thus  much  more  complicated  in  Eul'erian 
than  in  Lagrangian  coordinates.  Although  the  Eulerian  system  is 
preferable  to  the  Lagrangian  from  purely  dynamical  and  kinematical 
points  of  view,  the  Lagrangian  method  must  be  considered  the  simpler 
from  a  thermodynamical  viewpoint;  the  Eulerian  is  the  more  practi¬ 
cal  method  for  discussion  of  problems  which  are  mainly  dynamic, 
e.g.,  complicated  motions  of  simple  fluids,  the  Lagrangian  for  the 
study  of  problems  which  are  mainly  thermodynamic,  e.g.,  simple  motions 
of  media  having  complex  thermodynamic  properties.  If  the  amount  of 
water  vapor  is  assumed  to  be  constant  in  time,  and  no  fluid  or  solid 
water  is  supposed  to  be  present,  the  equation  of  state  for  the  air 
can  be  written 

ps  =  R@  .  (C-20) 


When  the  air  is  dry,  or  of  constant  composition,  the  individually 
constant  R  is  the  same  for  all  particles,  i.e.,  constant  in  time 
and  space.  In  this  case  C-20  can  be  used  directly  in  the  Eulerian 
system.  When  the  parameter  R  varies  from  one  particle  (a,  b,  c)  to 
another,  or  from  one  geometric  point  (x,  y,  z)  to  another,  R  is 
easily  eliminated  by  an  individual  logarithmic  differentiation  of 
Equation  C-20.  Hence, 


£ 

P 


+ 


(C-21) 


represents  a  more  general  form  of  the  atmospheric  equation  of  state. 


THE  ENERGY  EQUATIONS 

According  to  the  first  principle  of  thermodynamics,  the  change 
in  the  total  energy  of  any  fluid  particle  is  equal  to  the  work  done 
upon  it  by  all  forces,  external  as  well  as  internal,  plus  the  energy 
added  to  it  in  the  form  of  heat.  In  meteorology  we  neglect  for  the 
most  part  the  forms  of  energy  due  to  electric  and  magnetic  forces, 
and  consider  only  the  two  dynamically  and  thermodynamically  most 
important  types,  namely  the  kinetic  and  the  internal  energies.  The 
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kinetic  energy,  equal  to  1/2 v  per  unit  mass,  is  a  measure  of  the  in¬ 
tensity  of  the  macroscopic  motions  of  the  air  particles.  The  internal 
energy,  denoted  by  e,  is  an  expression  of  the  intensity  of  the  molecu¬ 
lar  motions ;  in  an  ideal  gas ,  this  energy  is  simply  proportional 
to  the  temperature.  The  heat  added  to  the  mass  unit,  by  processes 
of  radiation,  condensation,  etc.,  is  denoted  by  6w;  it  will  be  con¬ 
sidered  a  well-defined  function  of  time  and  space,  and  its  value 
per  unit  time  will  be  denoted  by  w  =  Dw/Dt.  The  principle  of  energy 
can  thus  be  written: 

(e  +  =  "  sv •  VP  -  p5  +  sv-S  +  S6  +  *  .  (C-22) 


This  equation  will  be  called  the  complete  energy  equation  or  the 
dynamical-thermodynamical  energy  equation,  since  it  contains  quan¬ 
tities  of  dynamical  as  well  as  thermodynamical  nature. 

A  special  dynamical  equation  of  energy,  independent  of  all 
transformations  of  heat  into  energy,  can  be  derived  from  the  equa¬ 
tion  of  motion  through  a  scalar  multiplication  by  v, 

(i?-2)  =  v?l  “  sv.vp  +  sv.S  .  (C-23) 


This  equation  simply  establishes  an  equality  between  changes  in  kinetic 
energy  and  the  external  work  done  by  gravity,  pressure,  and  friction, 
and  contains  neither  the  internal  work  of  pressure  and  friction,  nor 
the  internal  energy,  nor  the  added  heat.  If,  in  particular,  the  geo- 
pctential  0  at  a  fixed  geometric  point  is  independent  of  time,  we 
find 


D 

T5t 


sv.  (- vp  +  S) 


(C-24) 


In  this  form  the  dynamical  energy  equation  states  that  the  sum  of 
kinetic  and  potential  energies  is  changed  only  by  the  external  work 
done  by  the  forces  of  pressure  and  friction. 

Subtracting  Equation  C-23  from  C-22  now  gives  the  thermodynamical 
equation  of  energy, 
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w  +  s6  =  e  +  ps 


(C-25) 


or  if  e  =  cv0  and  the  potential  temperature  = 


9  are  introduced 


w  +  s6  = 


c  £  o 
P  <9 


(C-26) 


This  equation  connects  the  heat  supply  with  the  internal  work  of 
pressure  and  friction  and  with  changes  in  the  internal  energy,  but 
contains  neither  the  kinetic  energy  nor  the  external  work. 

Equations  C-22,  C-23,  and  C-24  which  establish  relations  be¬ 
tween  different  forms  of  energy,  do  not  say  that  one  definite  form 
of  energy  is  transformed  into  another  definite  form.  The  relations 
in  no  way  enable  us  to  follow  the  transformation  of  some  well-defined 
"amount11  of  energy.  Such  transformations  are  different  in  different 
motions,  and  a  complete  determination  for  a  given  motion  can  be 
carried  through  only  by  a  complete  knowledge  of  the  dynamics  of  the 
motion . 

If  the  energy  of  the  suspended  water  and  ice  is  negelected,  the 
internal  energy  of  moist  air  is  equal  to  cy9.  Since  the  quantity 
cv  varies  only  slightly  from  one  particle  to  another,  we  can  in 
practice  consider  it  constant  and  equal  to  the  specific  heat  cyd 
at  constant  volume  of  dry  air. 

The  heat  w  given  to  or  taken  from  the  air  particles  is  due  to 
process  of  evaporation,  condensation,  and  radiation,  and  also  to 
eddy  conductivity  if  we  consider  smoothed  motions.  If  we  assume 
the  composition  of  a  particle  to  remain  constant,  there  will  be  no 
heat  exchange  associated  with  evaporation  and  condensation.  There 
still  remains,  however,  the  radiative  heat  which  is  always  a  rather 
complex  quantity,  depending,  among  other  things,  on  the  absorptive 
power  of  the  air,  i.e.,  on  the  amounts  of  atmospheric  water  in  dif¬ 
ferent  phases.  It  seems  possible  to  discuss  the  large-scale  atmos¬ 
pheric  motions  by  means  of  the  several  equations  governing  this  type 
of  motion.  These  equations  form  a  complete  system  for  the  variables, 

including  the  three  velocity  components  v  ,  v  ,  v  ,  the  pressure  p, 

x  y  z 
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the  density  p  (or  the  specific  volume  s),  and  the  temperature  0. 
solution  of  this  system  will  contain  a  certain  number  of  arbitrary 
functions  and  constants.  Some  of  them  can  be  determined  by  the 
boundary  conditions,  whereas  others  are  given  by  the  initial  con¬ 
ditions  , 

v(x,  y,  z,  tQ)  =  vq(x,  y,  z)  , 

p(x,  y,  z,  tQ)  =  PQ(x,  y,  z)  , 

s (x ,  y,  z,  to)  =  sQ(x,  y,  z)  , 

e(x,  y,  z,  tQ)  =  e0(x,  y,  z)  . 

Here  tQ  is  a  certain  initial  time,  and  vq,  pQ,  sq,  and  0q  are  given 
functions,  to  be  determined  empirically  by  observations. 


ANNEX  D 


CONSERVATION  EQUATIONS  IN  PERTURBATION  FORM 

The  solution  of  the  nonlinear  hydrodynunical  equations  is  ex¬ 
tremely  difficult.  In  most  cases  it  is  therefore  necessary  to  sub¬ 
stitute  for  the  exact  system  of  equations  other  approximate  systems. 
The  method  of  perturbation,  introduced  systematically  in  hydro¬ 
dynamics  by  V.  Bjerknes  (1927),  enables  us  to  discuss,  by  means  of 
linear  differential  equations,  all  motions  that  lie  in  the  vicinity 
of  any  known  so-called  fundamental  motion.  The  difference  between 
the  fundamental  motion  and  the  motion  differing  only  slightly  from 
the  fundamental  motion  is  termed  the  motion  of  perturbation.  The 
latter  motion,  if  sufficiently  small,  satisfies  a  particularly 
simple  system  of  equations.  It  may  be  that  treatment  by  the  method 
of  perturbation  is  the  best  way  as  a  starting  point  of  dealing  with 
the  nonlinearity  of  advection  in  the  Eulerian  system.  The  following 
treatment  of  the  method  of  perturbation  as  applied  to  advection  of 
meteorologic  phenomena  is  that  given  by  Godske  et  al.  (1957),  who 
omitted  description  of  external  forces  (F^,  £em,  F^q11) ,  production 
and  loss  mechanisms,  and  heating  terms.  Inclusion  of  those  may, 
however,  be  treated. 

FUNDAMENTAL  MOTION  AND  PERTURBATION 

Let  the  velocity,  specific  volume,  pressure,  and  temperature 
characterizing  such  a  fundamental  motion  be  represented,  in  the 
Eulerian  system  of  coordinates,  by  the  quantities 


v  =  ±vx  +jvy  +  kvz,  s,  p,  e 


(D-l) 
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supposed  to  be  known  functions  of  the  independent  variables.  The 
velocity,  specific  volume,  pressure,  and  temperature  characterizing 
a  motion  that  differs  only  slightly  from  the  fundamental  one  can 
be  represented  by  the  quantities 

V  +  v  =  i(Vx  +  vx)  +  j(Vy  +  vy)  +  k(Vz  +vz)  , 

s  +  s?  p  +  p,  e  +  e  .  (d-2) 


Here  the  small  letters 

v  =  ivx  +  jvy  +  kvz,  s,  p,  e  (D-3) 

characterize  the  small,  unknown  "motion  of  perturbation" ;  when  added 
to  the  fundamental  motion  (Eq.  D-l) ,  they  are  determined  by  the 
above-discussed  hydrodynamical  equations.  If  the  perturbation  is 
small,  squares  and  products  of  the  quantities  D-3  are  supposed  to 
be  negligible,  so  that  a  "linear  system  of  equations  of  perturba¬ 
tion”  can  be  derived.  If  the  quantities  D-3  represent  a  solution 
of  this  system  approximate  solutions  of  the  hydrodynamical  equations 
are  thus  obtained  by  adding  D-l  and  D-3.  The  accuracy  of  the  approxi 
mation  improves  as  the  intensity  of  the  motion  of  perturbation  de¬ 
creases  . 

THE  GENERAL  EQUATIONS  OF  PERTURBATION 

Let  us  now  consider  the  equations  of  perturbation  for  a  friction 
less  fluid.  The  resultant  motion  must  satisfy  the  Eulerian  equations 
of  motion,  continuity,  state  and  thermodynamical  energy.  So  that 

~  +  +  y)  +  2fi  x  (y  +  v)  +  (S  +  s)  v(P  +  p) 

+  v($+0t)  =  O  (D-4) 

-  v-(v  +  v)  =  0  (D-5) 


S  +  s 


(S  +  s)  +  (y  +  y)-v(S  +  s) 
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x(s  +  s,  p  +  p,  0  +  e,  §  +  s,  ...)  =  o 


(D-6) 


W+ti=£+§+( P+  p)(£  +  S)  (D-7) 

where  the  perturbation  in  f  (produced,  tor  instance,  by  tidal  forces) 
is  denoted  by  0t  so  as  to  avoid  confusion  with  the  latitude. 

When  the  perturbation  is  produced  by  some  impulse  changing 
neither  the  external  forces  nor  the  heat  supply  to  a  physical  particle, 
the  quantities  V0j_  and  w  vanish;  such  a  perturbation  will  be  called 
"free."  If  the  perturbation  is  caused  by  a  small  additional  force 
whereas  the  heat  supply  to  the  particles  remains  unchanged,  we  have 
V0t  ^  0,  w  =  0.  Finally,  if  the  perturbation  is  produced  by  a  small 
additional  heat  supply  (e.g.,  a  daily  variation  of  insolation), 
we  have  w  ^  0,  v0t  =  0. 

Equations  D-4  to  D-7  inclusive  must  be  supplemented  by  boundary 
conditions.  Since  the  surface  of  discontinuity,  assumed  to  consist 
always  of  the  same  particles,  is  generally  also  disturbed,  its  equa¬ 
tion  in  the  motion  D-2  can  be  written 

F(x,  y,  z,  t)  +  f(x,  y,  z,  t)  =  0  (D-8) 

where  f  characterizes  the  perturbation  in  the  position  of  the  sur¬ 
face.  At  this  surface  the  kinematic  boundary  conditions 

f-?-  +  (V  +  y)-v(F  +  f)  =  0, 

a(F  +12  +  (V7  +  v')-V(F  +  f)  =  o  (D-9) 

must  be  satisfied,  as  well  as  the  dynamic  boundary  condition 

P  +  p-  P'-p'  =  0  (D-10) 

and  the  mixed  boundary  conditions 


65 


d(P  +  P  -  P —  P.  2  +  (y  +  ^).  v(P  +  p-  P'-p')=0 

at 

d(P  +  p  -  ?[  -  pJ2  +  (y '  +  ^')-v(P  +  p-  P'  -p')=0  .  (P-11) 

at  "" 

The  known  fundamental  motion  D-l  satisfies  a  corresponding 
system  of  equations,  obtained  simply  by  neglecting  all  small  letters 
in  equations  D-4  to  D-ll.  Thus  the  differential  equations  of  motion, 
continuity,  state,  and  energy  characterizing  the  fundamental  state 
become 


||  +  y.  vV  +  2&xV  +  S 


VP  +  y§  =  0 


||  +  V- VS  -  S  v-V  =  0 
at 

X(S,  P,  e,  s,  ...)  =  0 


w  -  £  -  pS  =  o  . 


The  undisturbed  surface  of  discontinuity  has  the  form 

F(x,  y,  z,  t)  =  0 


(D-12) 

(D-13) 

(D-14) 

(D-15) 

(D-16) 


and  the  corresponding  kinematic,  dynamic,  and  mixed  boundary  con¬ 
ditions  can  be  written 

||  +  V-vF  =  0,  ||  +  V'-vF  =  0  (D-17) 

©t  dt 

P  -  P'  =  0  (D-18) 

p/ -  +  v-v(P  -  p')  =  o 


a(P  -  p7) 

at 


+  v'-y(P  -  P')  =  o 


(P-19) 
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■s* 


The  system  of  Equations  D-12  to  D-15  can  now  be  subtracted 
from  the  system  D-4  to  D-7  valid  for  the  overall  resultant  motion. 
If  after  subtraction  squares  #nd  products  of  the  small  quantities 
characterizing  the  perturbation  are  neglected  (or  if  a  variation 
is  performed  on  the  equations  of  the  fundamental  motion) ,  the  fol- 
-otfing  linear  equations  of  perturbation  (linear  and  homogeneous  in 
v,  s,  p,  ...)  are  obtained: 


3y 

at 


+  V-  w+  y.yV  +  2QXv+S  vp  +  s  vP+  V0t  =  0 


(D-20) 


—■  +  v*ys  +  v-vs  -  s  v-y  -  s  v-y  =  o 

O  C 


Is  3  p  +le  (ft  + 


+  •  •  •  —  0 


(D-21) 


(D-22) 


I?  -  ft  '  P  ft  *  P  ft  +  ^<W  -  W  -  P®  -  Prf3) 


+  v-(sW  -  VE  -  PvS)  =  0  .  (D-23) 

•  • 

In  Equation  D-23  we  have  expanded  W  and  w,  although  they  are  not 
functions  of  state,  in  a  local  and  a  convective  part. 

The  equation  of  the  surface  of  discontinuity  is  of  the  form 

D-8 

F(x,  y,  z,  t)  +  f (x ,  y,  z,  t)  =  0  (D-24) 


thus  linear  in  f,  but  inhomogeneous.  The  formulas  representing  the 
kinematic  and  mixed  boundary  conditions,  however,  are  of  the  homo¬ 
geneous  type,  since  by  means  of  D-17  to  D-19  Equations  D-9  to  D-ll 
can  be  reduced  to 

M  +  v-vf  +  v-vF  =  0,  +  V'-vf  +  v'-yF  =  0  (D-25) 
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p  +  p  -  P' 


p  '  -  0 


(D-26) 


a(p  -  pO 

at 


+  y-7(p  -  p')  +  v*  v(P  -  P')  =  o 


d(pai~~  ^  +  -  p')  +  v'-vCP  -  P')  =  O  .  (D-27) 

The  system  of  differential  Equations  D-20  to  D-23  contains  the  po¬ 
tentials  $  and  0t,  assumed  to  be  known  field  functions,  constant  in 
time  or  variable,  and  the  internal  energies  E  and  e,  dependent  on 
the  temperature  0  and  its  perturbation  0.  Finally,  the  equations 
contain  the  amount  of  heat  W  given  to  the  particle,  as  well  as  its 
variation  w.  These  quantities  generally  depend  on  external  agencies; 
only  when  they  are  assumed  to  be  known  beforehand  will  the  system 
of  Equations  D-20  to  D-23  be  complete,  containing  just  as  many 
variables  as  equations. 

PERTURBATION  EQUATIONS  IN  SPHERICAL  COORDINATES 
Introducing 

V  V 

_  jd_  +  y  _9_  , _ \  5  _cg  _S_ 

Bt  r  ar  r  cos  co  r  3tp  ’ 

(at)s  -vr  a?  +  r  cos  cp  +  0^  (D-28) 

we  find  for  the  fundamental  state,  D-12  and  D-13,  the  following  equa¬ 
tions  of  motion  and  of  continuity  in  spherical  coordinates  in  the 
case  of  absolute  motion: 
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s  9vx  s  a 

r  cos  cp  J\  ~  r  cos  cp  (cos  'P  ^ 


ftV, 


r  cos  cp  ax 


r  cos  cp  5c p  ■  cp 


fr„  <cos  9  V  =  0 

(D-32) 


The  equations  of  state,  energy,  and  so  forth  are  easily  written  down 
by  means  of  the  general  expressions  for  the  fundamental  state  D-14 
to  D-19  and  for  the  perturbation  D-22  to  D-27. 
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ANNEX  E 


DIFFERENCE  EQUATIONS 


As  a  means  of  expressing  in  terms  of  discrete  steps  the  con¬ 
tinuous  functions  of  differential  equations,  we  resort  to  the  anal¬ 
ogous  difference  equations.  The  difference  equation  analogs  of 
differential  equations  may  take  a  variety  of  forms,  each  having 
different  ease  of  computation,  truncation  error,  and  stability  of 
solution  when  the  time  and  space  steps  are  made  to  be  small. 

If  the  discrete  space  positions  taken  by  the  space  steps  are 
denoted  by  subscripts  (i,  j,  k)  in  the  (x,  y,  z)  directions  and 
the  time  position  taken  by  the  time  step  is  denoted  by  the  super¬ 
script  (n)  in  the  time  dimension,  then  we  may  write  the  difference 
equation  analog  of  a  differential  by  reference  to  the  position  or 
time  for  which  the  evaluation  is  made.  Thus,  a  space  differential 
may  be  expressed  as  "forward,"  "centered,"  or  "backward,"  as  follows: 


Forward  space  difference  (in  x-dimension) 

(E-l) 


n  xi 

0  i+1/2  "  0  i-1/2 
AX 


(60)^ 


Centered  space  difference 
(in  x-dimension) 


(E-2) 


n 

0  i 


A-l 


Ax 


Backward  space  difference  (in  x-dimension) 


Similar  analogs  may  be  written  for  the  2  other  space  dimensions  and 
the  time  dimension. 


We  may  write  difference  equation  analogs  for  the  simple  dif¬ 
fusion  equation 


au 

St 


o  =  constant  >  0 


(E-4) 


as  an  explicit  difference  equation 


(a)  in  terms  of  a  "forward  time  difference,  centered  space 
difference"  form: 


n+1  n 

-3  — =  a  (62u)n.  =  a 

j 


,  n  n  n  n 

Uj+1  -2U  j  +  Uj-1\ 

(ax)2 


(E-5) 


(b)  or  a  "backward  time  difference,  centered  space  difference" 
form: 

un+1  -  ‘*n-  ,  nJ.,  un+J  -  2un+1  -  u"+* 

— - 1  =  o  (62u)^+1  =  a  -itL-J, - 2zl  {E.6) 


at 


(ax) 


or  as  a  combination  of  the  "backward"  and  "forward"  forms,  which  is 
an  implicit  scheme  for  evaluating  Eq.  E-4,  is  given  by 


un+1  - 

_3 _ 

at 


CT 


'  ,2  .n+1  .,2  xn 

0  (6  u)j  +  (1  -  0)(6  u)j 

(ax)2 


oseM. 

(E-7) 


Richtmyer  and  Morton  (1967)  have  described  a  review  of  tne 
several  approximations,  as  given  in  Table  E-l.  In  this  summary  of 
various  difference  systems  for  the  simple  diffusion  equation,  the 
diagram  at  the  left  indicates  the  type  of  differencing.  The  points 
shown  are  those  points  of  the  net  in  the  x,  t  plane  used  in  one 
application  of  the  formula;  the  t-axis  points  to  the  top  of  the  page 
and  x-axis  horizontally  to  the  right.  Where  three  points  are  con¬ 
nected  by  a  horizontal  line,  they  are  used  in  forming  the  second 
space  difference  ($  u) ,  and  where  two  points  are  connected  by  a 
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TABLE  E-l.  FINITE  DIFFERENCE  APPROXIMATIONS  TO 


(Special) 


»  |  *y'  3*  )*  ■ 

_J — 


'-f 


=  q.£-H  ,  a  =  constant  >  0 
at  » 

(after  Richtmyer  &  Morton  1967) 


un*l 


(HOj 

J(A*>* 
e  -  0(At)  ♦  0[<ax)2] 

explicit,  stable  if  oit/(A*)2  -  constant  t  1/2 


un*l 


-i  = « 


2(AX>* 

\ 

Crank  and  Nicholson  (1947) 
e  =  0[ (At)2}  +  0[(a*)2] 
implicit,  always  stable. 
n*l  ..n  ,  ,J..*n+l 


_J_  =  < 


<a‘U)3' 

(a*) 

Laasoner  (1949) 
e  =  0(aO  *  0[(AX)2] 
implicit,  always  stable. 

a.  Same  as  1,  but  with  oftt/(A x)‘  =  l/o 
•  Ot(lC)?]  =  OlCt.)4] 
special  case  of  1,  stable. 

S.  u"*1  -  u,  #<*2u)"*X  *  (l-«)(6‘a)n 


*'  l A* ) 

where  9  *  constant,  0  *  e  *  1 

e  =  0( At )  ♦  Or(Ak)2} 

for  0  s  e  <•  1/2,  stable  if  cbt /{&*)*  = 

■  or.start  *  1/(2  -  4a);  for  1/2  5  9  ?  1 
ilway".  stable 

inr  ljdes  i,  ....  4  as  jpecial  -*ses. 

5-<r«>  ss  but  with  4  =  1/7  -  ( A*  )2/l/cftt 
*>  -  7'  (At  >2]  ♦  ( A*  )*  ] 

'•table . 


- Tit 

always  wstable. 


*>♦1  r—  1 


~7I7" 


(•x> 

•here  .»/a*  - 

1  ( At  )*' j  ♦  OSAO's  *  7  Kir 


3/2 

_  1  /- 


7  At  7  At  (ax)' 


e  =  0[(At)  ]  ♦  O'(ftx)  ) 

always  stable. 

(in*l  tin  un  _  un-l 


“•  11  *  ” — 3—  ■  — w 

(A*)? 

where  9  =  constant  t  0.  fJAt/fftx)* 
e  =  0(/»t )  ♦  0;  (ax)*] 
always  stable 

contains  5,  9  as  special  rases. 
11.  Saxe  as  10  but  with 


3/2  -(A*)  /1206t 

-2/2*(flx)2/22o£t 


'  hx; 


a  _  1  4  (A*)' 

7  12ofit 

,  =  O'(lt)2'.  =  0[(,K)*1 
always  stable. 

n*l  r.  ,n+l 

1  u .  . ,  -  u .  ,  r  'J 


17  ■ 


"'ft!  ” 


n+1  n 

i  11  i- 1  '  :-i 

TT  - ^  ' 


e  ■-  0’(aI>*]  ♦  0r(A>.)  * 
always  stable 
Identical  with  6. 


23.  t  2/2  ^  -  2^.2  ♦ 


1  1  1 
IT  5  TT 


VJU^1  -  .->0°  ♦  1/.  /-1 


1  ^f::i  -  -2-1  *  i/“ 

’  TT  - ‘ - It - 


(  6*  a  ), 


n  OK  At)*  ;  *  0[(A»  )  ; 
always  .tulle. 


/n  *i  n  *i  ,n*l  .,r*l» 

!VV'*1  ’  ’  **:  ’ 


/(.*«>' 


Su  Ill' 


lu  For* 
»pii  »' 


F ranfre  j  fir.') 
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vertical  line  they  are  used  in  forming  the  time  difference.  When 
there  are  two  or  more  sets  of  points  of  either  kind,  the  weights  used 
are  indicated  at  the  side  or  underneath.  The  table  also  shows  the 
order  of  the  truncation  error  (e)  and  the  conditions  for  stability  of 
solution. 

The  5-point  explicit  scheme  (#1  of  Table  E-l)  with  forward  time 
difference  is  that  employed  by  Leith  (1965)  and  to  be  used  for  the 
S-M-T  Model. 

Richtmyer  and  Morton  (1967)  point  out  that  in  general  the 
stability  of  an  approximation's  solution  is  determined  by  the  sta¬ 
bility  of  the  solution  of  the  highest  order  derivative. 

As  a  means  of  replacing  a  complicated  (say  multidimensional) 
problem  by  a  succession  of  simpler  (say  1-dimensional)  ones,  the 
methods  of  splitting,  or  methods  of  fractional  steps,  were  devised 
by  Soviet  mathematicians  (Bagrinovski  and  Godunov  1957,  B'yakcnov, 
1962,  Marchuk  1964).  These  have  been  applied  by  Leith  (1965)  and 
are  proposed  for  the  S-M-T  Model  in  the  following  form: 


,n+i/2  = 


...n 


) 2  k 


2yn  .  +  f 

J.K 


(E-8) 


n+1 

V 


,n+l/2 

iik 


n+1 

j+l.k 

J  i 


n+1/  2 
j+l-jk 


2wn+1/2  +  yn+1/2\ 

(E-9) 


The  computational  process  involved  in  a  single  time  step  is  divided 
into  two  cycles.  In  the  first  of  these,  advection  in  the  y-direction 
only  is  evaluated;  whereas  in  the  second,  advection  in  the  x-direction 
only  is  evaluated  starting  of  course  with  the  results  of  the  first 
cycle.  The  fractional  time  step  notation  is  a  convenience,  but  the 
results  of  the  first  cycle  have  no  particular  physical  significance. 
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CONSIDERATIONS  FOR  CHOICE  OF  SPACE  AND  TIME  STEPS 


ALTITUDE,  km 


PER  (hr)  100  10 

♦ 

4d 


0.1  0.01  0.001 

♦ 

1  min 


4.0 

3.0 

<N 

I 

O 

S  2.0 

N 

E 

1.0 

0 

lyr  l/2yr  30d  lOd  6d  3d  24hr  12hr  6hr  3hr  2hr  Ihr  20m  10m  5m  2m  !m 

SPECTRUM  WIND  SPEED  AT  CARIBOU,  MAINE  (After  Oort  and  Taylor,  1969) 


FIGURE  F-4.  Horizontal  Wind  Speed  Spectra 
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STI 1-20-68-3  PRESSURE  (mb)  OR  DENSITY  (9  m"3) 


FIGURE  F-5.  Temperature,  Pressure,  Density,  Molecular  Weight  (  U.  S. 
Standard  Atmosphere,  1962 )( After  Valley,  1965) 
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ANNEX  G  ; 

HEAT  SOURCES  AND  SINKS 

•  Domdin 

Heat  Source/Sink 

0-20 

Emission  from  surface 

Emission  by  H20  vapor  (50^) 

Emission  by  C02  (15^,) 

Absorption  by  H20  vapor 

Convection/Conduction  upward 

20-50 

Emission  by  C02  (IR) 

Emission  by  H20  i 

Absorption  by  0^  (  ) 

Conduction  (down) 

Absorption  by  02>  H20,  N^O ,  CH4 

50-30 

Emission  by  C02  ( IS p,) 

Emission  by  0^  (9.6g,) 

Emission  by  0  (63p,) 

Absorption  by  0^  (uv) 

Absorption  by  02  (uv) 

Particle  precipitation 

Conduction  (upward) 

80-150 

Emission  by  0  (63^) 

Absorption  by  0^  (uv) 

Absorption  by  C>2,  N2  (xr) 

Particle  precipitation 

Conduction  (down) 

Joule  heading  by  ionosphere  dynamics 

Dissipation  of  gravity,  tidal,  Rossby  waves 

Recombination  of  down-t ransported  0 

150-300 

Absorption  by  N9,  0?  (Euv) 

Particle  precipitation 

_ 

± 


ANNEX  H 


i 


CHEMICAL  REACTIONS  FOR  PRODUCTION  AND  LOSS 
OF  ATMOSPHERIC  CONSTITUENTS 

The  principal  production  and  loss  mechanisms,  with  greatest 
time  and  space  variability,  are  those  due  tp  ionization  and  dis¬ 
sociation  resulting  from  the  solar  electromagnetic  radiation,  which 
is  wavelength  and  flux  dependent. 

The  dissociation  energies  and  ionization  potentials  of  atmos¬ 
pheric  species  of  interest,  as  listed  by  Bortner  and  Kummler  (1968), 
are  given  in  Tables  H-l,  H-2,  and  H-3.  Since  the  more  intense  X-ray 
radiation  may  produce  several  icn-pairs  per  quantum,  the  more  ener¬ 
getic  radiation  is  expected  to  produce  ion-pairs  at  a  rate  propor¬ 
tional  to  the  next  smaller  integer  value  of  ( J/IP) ,  wherein  J  is 
-1  3 

the  ergs  sec  deposited  in  each  cm  of  the  radiation  path,  and  IP 
is  the  ionization  potential  of  the  species. 

Other  important  sources  of  production  and  loss  of  chemical  con¬ 
stituents  are  the  chemical  processes,  defined  by  Bortner  and  Kummler 
(1968),  which  are  described  in  Figs.  H-l  to  H-31  inclusive.  Only 
the  most  important  reactions  and  their  rate  coefficients  are  shown. 

The  species  shown  in  the  Figs.  H-l  to  H-31  are  in  the  following 
order: 

a.  Major  neutral  species:  H^O,  CO^,  0^  (Fig.  5-21) 

b.  Minor  neutral  species:  0,  0(1D),  02(1A),  O^1^),  0^,  N, 

no,  no2,  n2o,  h,  h2,  oh,  ho2,  h2o2 

c.  Positive  ions:  N*,  N+,  0+,  0*,  N0+ 

d.  Negatively  charged  species:  e,  02,  0  ,  0^,  00^,  N02,  N0^ 

Preceding  page  blank 
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e.  Metallic  species: 
Mg  and  Cu. 


of  form  M,  MO,  M02> 


M  of  metals  Na, 


For  each  of  the  reactions  shown,  a  rate  coefficient  is  given 

2  -1 

in  cgs  units,  that  is  in  cm  sec  per  particle  for  two  body  and  in 
6  -1  2 

cm  sec”  per  (particle)  for  three  body  reactions.  Particle  den¬ 
sities  are  in  cm-^  and  temperature  in  °K.  M  is  used  to  represent  a 
collision  partner.  All  values  for  rate  constants  are  given  by  the 
constants  a,  b,  and  c  for  a  rate  constant  of  the  form 

\b 


R  =  a 


(  T  V 

\3UUJ 


exp  (-c/T). 


The  rate  coefficients  may  have  error  of  factors  2  to  10. 


Bortner  and  Kummler  (1968)  list  33  atmospheric  species,  with  ap¬ 
proximately  8  reactions  each,  or  244  reactions. 


TABLE  H-l. 

DISSOCIATION  ENERGIES, 

IONIZATION 

POTENTIALS,  AND 

ELECTRON  AFFINITIES  OF 

SPECIES  OF 

INTEREST 

D(ev) 

EP(ev) 

EA(ev) 

O  ” 

- 

13.618 

1.478 

N 

- 

14.532 

o„ 

5.115  (O-O) 

12.063 

0.43 

N2 

9,759  (N-N) 

15.580 

NO 

6.507  (N-O) 

9.267 

0.3 

03 

1.051  (0-02) 

12.80 

1.9 

no2 

3. 11G  (O-NO) 

9.78 

4.0 

n2o 

1.677  (N2-0) 

12.894 

nc>3 

2.122  (O-NO2) 

Ar 

- 

15.759 

H 

- 

13.598 

0.754 

H2 

4.477  (H-H) 

15.425 

OH 

4.395  (0-H) 

13.34 

1.83 

H„0 

5.116  (H-OH) 

12.619 

H&2 

2.04  (H-02) 

H2Ojj 

2.22  (HO-OH) 

CO 

11.108  (C-O) 

14.013 

C02 

5.453  (O-CO) 

13.769 

C03 

Li 

5.390 

Na 

5.138 

K 

4.339 

Cs 

3.893 

Mg 

7.644 

Ca 

6.111 

8a 

5.210 

Si 

8.149 

be 

7.87 

Bortner  and  Kummler:  "The  Chemical  Kinetics  and  the  Composition  of 

the  Earth’s  Atmosphere,"  G.E.  Reentry  Systems 
Science  Report  GE-9500-ECS-SP-1,  24  July  1968. 
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TABLE  H-2.  WAVELENGTHS  CORRESPONDING  TO  IONIZATION  POTENTIALS,  A 


N2 

796 

N 

853 

0 

911 

H 

912 

°2 

1029 

no2 

1269 

NO 

1339 

Si 

1547 

Fo 

1601 

Mg 

1649 

Ca 

2062 

Na 

2453 

K 

2905 

TABLE  H-3.  WAVELENGTHS  CORRESPONDING  TO  DISSOCIATION  ENERGIES 


D(ev) 

VA> 

o,  - 

0(3P)  +  0(3P) 

5.115 

2425 

°22'- 

OCP)  +  0(1D) 

7.082 

1751 

X 

NO  - 

N(4S)  +  0(3P) 

6.507 

1907 

NO  - 

N(4S)  +  O^D) 

8.474 

1464 

NO  - 

N(2D)  +  0(3P) 

8.891 

1395 

X 

°q  - 

o(3p)  +  o2(2r) 

1.051 

11804 

03  - 

O(Jd)  +  o2(jE) 

3.018 

4111 

°3  - 

Of3?)  +  Og^A) 

2.028 

6117 

O 

°3  * 

O(Jd)  +  o2(Ja) 

3.995 

3105 

X 

O 

°3  " 

om  +  o2{1£) 

2. 678 

4633 

°3  - 

0(XD)  +02(1^) 

4.645 

2671 

X 

no2- 

Of1D)  +  NO^it) 

3.116 

3981 

X 

Kffocllvo 

1  diaaociation  procosaoa 
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FIGURE  H-l .  Chemical  Kinetics  of  Otygen  Molecules 
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FIGURE  H-4.  Chemical  Kinetics  of  Oxygen  Atoms 


Chemical  Kinetics  of  Exc 


+  h  v 


FIGURE  H-7.  Chemical  Kinetics  of  (b  Lg) 
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FIGURE  H-8 


FIGURE  H-9.  Chemica 
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FIGURE  H-10.  Chemical  K 
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FIGURE  H-l 


FIGURE  H-13.  Chemical  Kinetics  of  Nitrous  Oxide 


FIGURE  H-15.  Chemical  Kinetics  of  Hydrogen  Molecules 


FIGURE  H-17.  Chemical  Kinetics  of  HO„  Radicals 


FIGURE  H-18.  Chemical  Kinetics  of  Hydrogen  Peroxide 


FIGURE  H-19.  Chemical  Kinetics  of  Molecular  Nitrogen  Ions 


FIGURE  H-20. 
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FIGURE  H-22.  Chemica 


FIGURE  H-23.  Chemical  Kinetics  of  Nitric  ^cide  Ions 
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FIGURE  H-24.  Chemical  Kinetics  of  Free  Electrons 


FIGURE  H-26.  Chemical  Kinetics  of  Negative  Atomic  Oxygen  Ions 


FIGURE  H-27.  Chemical  Kinetics  of  Negative  Ozone  Ions 


FIGURE  H-28.  Chemica 
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FIGURE  H-29.  Chemical  Kinetics  of  Negative  Nitrogen  Dioxide  Ions 


FIGURE  H-30.  Chemical  Kinetics  of  Negative  Nitrogen  Trioxide  Ions 
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FIGURE  H-31.  Chemical  Kinetics  of  Metal  Atoms,  Metal  Oxide 
Molecules  and  Metal  Ions  of  Na,  Mg,  Cu. 


ANNEX  I 


RADIATION  ENERGY  AND  HEAT  EQUATIONS 

INTRODUCTION 

The  Stratosphere-Mesosphere-Thermosphere  as  given  in  Appendix  E 
is  based  on  the  grid  structure  and  diff erenceing  schemes  of  the 
Leith  Tropospheric  model,  with  radiation  energy,  heat  equations, 
electromagnetic  effects ,  etc . ,  as  added  features .  In  order  to  get 
some  idea  of  the  extra  computational  load  imposed  by  the  added 
features,  each  was  to  be  assessed  in  order.  After  looking  at  the 
radiation  energy  and  heat  equations,  it  became  apparent  that  the 
remaining  effects  would  not  contribute  significantly  to  the  load 
and  their  detailed  analysis  was  abandoned. 

This  Annex  presents  the  analysis  of  the  radiation  energy  and 
heat  equations. 

EQUATIONS 

The  radiation  energy  equation  uses  as  its  variable  input  data 

the  values  I  ,  which  are  received  from  the  SOLRAD  HI  satellite, 
m 

These  values  are  assumed  to  change  hourly.  The  equation  convolves 
this  information  along  with  other  known  information  about  the  con¬ 
stituents  of  the  upper  atmosphere  to  produce  coefficients  for  a  set 
of  simultaneous  linear  equations.  The  solution  of  these  equations 
produces  the  number  density  of  each  constituent  of  the  atmosphere 
at  each  designated  level  in  the  model. 

The  stratosphere  is  divided  into  6  layers,  the  mesosphere  into 
8  layers,  the  lower  thermosphere  into  6  layers,  and  the  upper  thermos¬ 
phere  into  8  layers,  for  a  total  of  28  layers. 

The  data  sets  involved  in  the  radiation  energy  equation  and 
their  origin  or  disposition  are  given  in  the  following  list: 

Preceding  page  blank  U9 


a-  -(or  a  - 

ui,3 

F(X)  - 

I  - 
m 

IP.  - 


SJ 


t,m 

S. 


n 


N.  - 
J 

G  - 


R.  - 


Given  array  (time  independent) 

Given  table  (30  to  100  values) 

SOLRAD-HI  input  values 

Given  constants 

Intermediate  results 

Intermediate  results  (save  for  heat  equation) 

Input  to  simultaneous  equations  (see  Eq.  1-1 
below) 

Results  of  simultaneous  equations 
Recursion  formula 

Matrix  of  values  selected  from  S  and  R 

l  l 

using  Annex  H 

Given  constants  that  go  into  G-matrix. 


The  equations  to  be  solved  and  the  definition  of  the  indices  are 
given  below  in  Table  1-1.  The  indices  are  as  follows: 


TABLE  1-1.  EQUATION  INDICES 


Domain 

Computation 

Cycle 

Limits 

t 

t-1 

Components 

3 

<:  42 

Wave  bands 

i 

s  20 

Altitude  levels 

k 

k  -  1 

31 

Equations,  based  on  Equations  32  and  33  of  this  Appendix  are  the 
following: 


S 


1 

IF" 

l 


J 


(1-1) 
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Solution  of  i  x  l 


n  =  simultaneous 
l 

equations 


(1-2) 


NjCk)  =  N..(h  -  1)  +  nj(k  -  1)  ,  (1-3) 

where  n  (k)  are  the  results  of  the  solution  of  the  i  %  l  simul- 
taneous  equation  set.  (See  Equation  37  of  this  Appendix.) 


Qt+1  =  Cf  +  At 


(1-4) 


TIMING  ANALYSIS 

The  indices  involved  are  given  more  precisely  in  Table  1-2: 


TABLE  1-2.  INDICES  BY  DOMAIN 


Number  of 
Components 


Wave  Bands 
Range  of  m 


Number  of 
Wave  Bands 


(j  and  /.) 

(nv,  to  M) 

(i  and 

Stratosphere 

18 

11  to  20 

10 

Mesosphere 

30 

11  to  20 

10 

Lower  Thermosphere 

42 

1  to  16 

16 

Upper  Thermosphere 

30 

8  to  11 

4 

Equation  1-1  may  be  broken  down  into  6  logical  computation  sequences: 

a  =  ZL  cj-:  -W*  “  D 


b  =  F(X)  •  a 
c  =  Exponential  of  b 


I  •  c 
m 
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The  total  computation  for  the  equation  can  be  determined  by 
substituting  the  proper  integer  indices  for  i,  m,  and  j  and  the 
number  of  microsec  for  each  of  the  6  computational  sequences  in  the 
following  formula: 

Total  Computation  =  -t[m(j  *a+b+c+d+c)+f]. 

However,  it  is  more  appropriate  to  determine  the  time  required  per 
point,  since  all  previous  estimates  and  timings  are  made  on  the 
basis  of  a  point.  This  posed  a  small  problem,  which  was  solved  by 
calculating  a  weighted  average  for  l,  m,  and  j,  determined  by  the 
number  of  layers  in  each  of  the  domains.  Thus  the  weighted  average 
for  l  and  j  was  30,  and  for  m  it  was  9. 

The  CDC-7600  was  used  as  the  basis  for  the  execution  time  per 
section,  because  of  a  good  timing  estimate  given  for  the  Leith 
troposphere  model  which  is  running  on  the  7600  at  the  Lawrence 
Radiation  Laboratory  in  Livermore,  California 

Without  overhead,  the  6  sections  were  timed  as  follows: 


Section 


Number  of 
CDC-7600 
Cycles 


a 

b 

c 

d 

e 

f 


19 

77 

95 

15 

23 

32 
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Using  these  values  and  the  weighted  averages  for  £,  m,  and  j ,  we 
have: 

Time  per  point  =  27.5  ns  per  cycle  x  30[9(30.19  +-  210)  +  32] 

=  5.8  ms  per  point. 

Equation  1-2  can  be  divided  into  2  logical  sections:  (a)  location 
of  the  coefficients  from  given  and  calculated  values  and  their  for¬ 
mation  into  a  matrix  suitable  for  solution,  and  (b)  solution  of  the 
simultaneous  equations. 

The  first  of  these  sections  consists  of  the  manipulation  of 
indexes  to  select  the  proper  coefficients.  A  set  of  predetermined 
pointers  cam  be  developed  for  each  domain  to  minimize  the  complexity 
of  this  process.  Without  showing  the  detail,  this  section  required 
approximately  0.15  ms  per  point. 

The  solution  of  the  simultaneous  equations  was  based  on  a  popu¬ 
lar  subroutine  used  at  LRL  on  the  CDC-6600.  For  actual  matrices  of 
equivalent  density  but  for  a  10  x  10,  the  time  required  on  the  6600 
was  approximately  1.0  millisec.  Since  the  solution  time  using  this 
routine  goes  up  approximately  linearly  as  a  function  of  the  number 
of  nonzero  values  in  the  matrix,  we  can  extrapolate  that  it  will 
require  about  9  millisec  for  the  30  x  30  matrix,  which  is  the  weighted 
average  size  for  the  S-M-T  Model.  Using  a  factor  of  5:1  for  the 
CDC-7600  over  the  CDC-6600  results  in  an  estimate  of  1.8  millisec 
for  solution  on  the  CDC-7600. 

Equation  1-3  is  simple  and  requires  an  average  of  31  iterations 
for  an  estimated  25  microsec  of  CDC-7600  computer  time. 

Equation  1-4,  the  heat  equation,  is  simple  by  virtue  of  the 
fact  that  partial  results  of  earlier  radiation  calculations  can  be 
used  to  achieve  this  result.  The  equation  is  estimated  to  require 
approximately  20  microsec  of  CDC-7600  time. 

The  total  time  for  Equations  1-1  through  1-4  is  found  to  be 
about  7.75  millisec  per  point  without  considering  overhead  and 
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bookkeeping  instructions-  Since  maximum  times  have  been  used  in 
all  of  the  previous  analyses,  it  is  felt  that  the  overhead  will  be 
more  than  offset  by  the  difference  that  would  occur  from  overlap 
of  arithmetic  instructions. 

The  Troposphere  hydrodynamic  model  after  Leith,  which  is  being 
run  on  the  CDC-7600  at  the  Lawrence  Radiation  Laboratory,  Livermore, 
has  been  estimated  to  take  about  180  microsec  per  point.  Thus,  from 
a  comparative  standpoint,  the  radiation  energy  and  heat  equations 
would  take  approximately  43  times  as  much  time  as  is  required  for 
the  hydrodynamical  computations. 

Since  it  was  determined  at  this  point  that  the  desired  result 
had  been  realized,  no  further  analyses  were  made  on  the  remaining 
features  of  the  S-M-T  Model.  They  were,  however,  estimated  to  be 
small  in  comparison  to  the  radiation  energy  calculation. 
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